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Preface 


The beauty of plants has attracted the attention of mathematicians for 
centuries. Conspicuous geometric features such as the bilateral sym- 
metry of leaves, the rotational symmetry of flowers, and the helical 
arrangements of scales in pine cones have been studied most exten- 
sively. This focus is reflected in a quotation from Weyl [159, page 3], 
“Beauty is bound up with symmetry.” 

This book explores two other factors that organize plant structures 
and therefore contribute to their beauty. The first is the elegance and 
relative simplicity of developmental algorithms, that is, the rules which 
describe plant development in time. The second is self-similarity, char- 
acterized by Mandelbrot [95, page 34] as follows: 


When each piece of a shape is geometrically similar to the 
whole, both the shape and the cascade that generate it are 
called self-similar. 


This corresponds with the biological phenomenon described by Herman, 
Lindenmayer and Rozenberg [61]: 


In many growth processes of living organisms, especially of 
plants, regularly repeated appearances of certain multicel- 
lular structures are readily noticeable..... In the case of a 
compound leaf, for instance, some of the lobes (or leaflets), 
which are parts of a leaf at an advanced stage, have the 
same shape as the whole leaf has at an earlier stage. 


Thus, self-similarity in plants is a result of developmental processes. 
By emphasizing the relationship between growth and form, this book 
follows a long tradition in biology. D’Arcy Thompson [143] traces its 
origins to the late seventeenth century, and comments: 


Organic form itself is found, mathematically speaking, to be 
a function of time.... We might call the form of an organism 
an event in space-time, and not merely a configuration in 
space. 


This concept is echoed by Hallé, Oldeman and Tomlinson [58): 


The idea of the form implicitly contains also the history of 
such a form. 
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The developmental processes are captured using the formalism of 
L-systems. They were introduced in 1968 by Lindenmayer [82] as a 
theoretical framework for studying the development of simple multicel- 
lular organisms, and subsequently applied to investigate higher plants 
and plant organs. After the incorporation of geometric features, plant 
models expressed using L-systems became detailed enough to allow the 
use of computer graphics for realistic visualization of plant structures 
and developmental processes. 

The emphasis on graphics has several motivations. A visual compar- 
ison of models with real structures is an important component of model 
validation. The display of parameters and processes not observable di- 
rectly in living organisms may assist in the analysis of their physiology, 
and thus present a valuable tool for educational purposes. From an 
aesthetic perspective, plants present a wealth of magnificent objects 
for image synthesis. The quest for photorealism challenges modeling 
and rendering algorithms, while a departure from realism may offer a 
fresh view of known structures. 

The application of computer graphics to biological structures is only 
one of many factors that contribute to the interdisciplinary character 
of this book. For example, the notion of L-systems is a part of formal 
language theory, rooted in the theory of algorithms. The application of 
L-systems to plant description has been studied by biologists, and in- 
volves various methods of general mathematics. Self-similarity relates 
plant structures to the geometry of fractals. Computer-aided visual- 
ization of these structures, and the processes that create them, joins 
science with art. 

The study of an area that combines so many disciplines is very stim- 
ulating. Some results may be of special interest to students of biology 
or computer graphics, but a much wider circle of readers, generally in- 
terested in science, may find mathematical plant models inspiring, and 
the open problems worth further thought. Consequently, all basic con- 
cepts are presented in a self-contained manner, assuming only general 
knowledge of mathematics at the junior college level. 

This book focuses on original research results obtained by the au- 
thors in the scope of the cooperation between the Theoretical Biology 
Group, directed by Aristid Lindenmayer at the University of Utrecht, 
and the Computer Graphics Group, working under the supervision of 
Przemyslaw Prusinkiewicz at the University of Regina. Technically, the 
book evolved from the SIGGRAPH ’88 and ’89 course notes Linden- 
mayer systems, fractals, and plants, published by Springer-Verlag in 
the series Lecture Notes in Biomathematics [112]. The present volume 
has been extended with edited versions of recent journal and conference 
papers (see Sources), as well as previously unpublished results. 

Aristid Lindenmayer is the author of the notion of L-systems which 
forms the main thread of the book. He also played an essential role 
in the reported research by suggesting topics for study, guiding the 
construction of specific plant models, monitoring their correctness and 
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participating in many discussions of biological and mathematical prob- 
lems. Seriously ill, Professor Lindenmayer co-authored and edited sev- 
eral chapters, but was not able to participate in the completion of this 
work. If any inaccuracies or mistakes remain, he could not prevent 
them. Still, in spite of unavoidable shortcomings, we hope that this 
book will convey his and our excitement of applying mathematics to 
explore the beauty of plants. 
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Chapter 1 


Graphical modeling using 
L-systems 


Lindenmayer systems — or L-systems for short — were conceived as 
a mathematical theory of plant development [82]. Originally, they did 
not include enough detail to allow for comprehensive modeling of higher 
plants. The emphasis was on plant topology, that is, the neighborhood 
relations between cells or larger plant modules. Their geometric aspects 
were beyond the scope of the theory. Subsequently, several geometric 
interpretations of L-systems were proposed with a view to turning them 
into a versatile tool for plant modeling. Throughout this book, an 
interpretation based on turtle geometry is used [109]. Basic notions 
related to L-system theory and their turtle interpretation are presented 
below. 


1.1 Rewriting systems 


The central concept of L-systems is that of rewriting. In general, rewrit- 
ing is a technique for defining complex objects by successively replacing 
parts of a simple initial object using a set of rewriting rules or produc- 
tions. The classic example of a graphical object defined in terms of 
rewriting rules is the snowflake curve (Figure 1.1), proposed in 1905 by 
von Koch [155]. Mandelbrot [95, page 39] restates this construction as 
follows: 


One begins with two shapes, an initiator and a generator. 
The latter is an oriented broken line made up of N equal 
sides of length r. Thus each stage of the construction begins 
with a broken line and consists in replacing each straight 
interval with a copy of the generator, reduced and displaced 
so as to have the same end points as those of the interval 
being replaced. 


Koch 
construction 


Grammars 


L-systems 


2 Chapter 1. Graphical modeling using L-systems 


initiator 


=_ oo 


generator 


Figure 1.1: Construction of the snowflake curve 


While the Koch construction recursively replaces open polygons, rewrit- 
ing systems that operate on other objects have also been investigated. 
For example, Wolfram [160, 161] studied patterns generated by rewrit- 
ing elements of rectangular arrays. A similar array-rewriting mecha- 
nism is the cornerstone of Conway’s popular game of life [49, 50]. An 
important body of research has been devoted to various graph-rewriting 
systems [14, 33, 34]. 

The most extensively studied and the best understood rewriting sys- 
tems operate on character strings. The first formal definition of such a 
system was given at the beginning of this century by Thue [128], but 
a wide interest in string rewriting was spawned in the late 1950s by 
Chomsky’s work on formal grammars [13]. He applied the concept of 
rewriting to describe the syntactic features of natural languages. A 
few years later Backus and Naur introduced a rewriting-based notation 
in order to provide a formal definition of the programming language 
ALGOL-60 [5, 103]. The equivalence of the Backus-Naur form (BNF) 
and the context-free class of Chomsky grammars was soon recognized 
[52], and a period of fascination with syntax, grammars and their appli- 
cation to computer science began. At the center of attention were sets 
of strings — called formal languages — and the methods for generating, 
recognizing and transforming them. 

In 1968 a biologist, Aristid Lindenmayer, introduced a new type of 
string-rewriting mechanism, subsequently termed L-systems [82]. The 
essential difference between Chomsky grammars and L-systems lies in 
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Regular 
Context-free 


Context-sensitive 





Recursively enumerable 


Figure 1.2: Relations between Chomsky classes of languages and language 
classes generated by L-systems. The symbols OL and IL denote language 
classes generated by context-free and context-sensitive L-systems, respec- 
tively. 


the method of applying productions. In Chomsky grammars produc- 
tions are applied sequentially, whereas in L-systems they are applied 
in parallel and simultaneously replace all letters in a given word. This 
difference reflects the biological motivation of L-systems. Productions 
are intended to capture cell divisions in multicellular organisms, where 
many divisions may occur at the same time. Parallel production ap- 
plication has an essential impact on the formal properties of rewriting 
systems. For example, there are languages which can be generated 
by context-free L-systems (called OL-systems) but not by context-free 
Chomsky grammars [62, 128] (Figure 1.2). 


1.2 DOL-systems 


This section presents the simplest class of L-systems, those which are 
deterministic and context-free, called DOL-systems. The discussion 
starts with an example that introduces the main idea in intuitive terms. 

Consider strings (words) built of two letters a and b, which may 
occur many times in a string. Each letter is associated with a rewriting 
rule. The rule a — ab means that the letter a is to be replaced by 
the string ab, and the rule b — a means that the letter b is to be 
replaced by a. The rewriting process starts from a distinguished string 
called the axiom. Assume that it consists of a single letter b. In the 
first derivation step (the first step of rewriting) the axiom b is replaced 


Example 


L-system 


Derivation 
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pee 


a 


b 
_! | 
aba 
tlh Ee 
abaab 


af 1b 


abaababa 


Figure 1.3: Example of a derivation in a DOL-system 


by a using production b — a. In the second step a is replaced by ab 
using production a — ab. The word ab consists of two letters, both of 
which are simultaneously replaced in the next derivation step. Thus, a 
is replaced by ab, 6 is replaced by a, and the string aba results. In a 
similar way, the string aba yields abaab which in turn yields abaababa, 
then abaababaabaab, and so on (Figure 1.3). 

Formal definitions describing DOL-systems and their operation are 
given below. For more details see [62, 127]. 

Let V denote an alphabet, V* the set of all words over V, and 
V+ the set of all nonempty words over V. A string OL-system is an 
ordered triplet G = (V,w,P) where V is the alphabet of the system, 
w € V* is a nonempty word called the aziom and P CV x V* isa 
finite set of productions. A production (a,y) € P is written as a > 
y. The letter a and the word y are called the predecessor and the 
successor of this production, respectively. It is assumed that for any 
letter a € V, there is at least one word y € V* such that a > yx. If 
no production is explicitly specified for a given predecessor a € V, the 
identity production a — ais assumed to belong to the set of productions 
P. An OL-system is deterministic (noted DOL-system) if and only if 
for each a € V there is exactly one y € V* such that a — y. 

Let = a,...@, be an arbitrary word over V. The word v = 
X1---Xm € V* is directly derived from (or generated by) jz, noted ps => 
v, if and only if a; — x; for alli =1,...,m. A word v is generated by 
G in a derivation of length n if there exists a developmental sequence of 
Words ig, fis... hy: Such that fig = Wy ly =v and (oS a Se => 


Ln- 
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Figure 1.4: Development of a filament (Anabaena catenula) simulated using 
a DOL-system 


The following example provides another illustration of the operation of 
DOL-systems. The formalism is used to simulate the development of a 
fragment of a multicellular filament such as that found in the blue-green 
bacteria Anabaena catenula and various algae [25, 84, 99]. The symbols 
a and b represent cytological states of the cells (their size and readiness 
to divide). The subscripts | and r indicate cell polarity, specifying the 
positions in which daughter cells of type a and b will be produced. The 
development is described by the following L-system: 


Wt Gp 
Pit Qp— aid, 
p2: a bia, (1.1) 


P3: bp — Qr 
pa: ba 


Starting from a single cell a, (the axiom), the following sequence of 
words is generated: 


iy 
a,b, 

bja,Q, 
ajajb,ajb, 

bj, A,b,0,A,-0, Ay Ay 


Anabaena 


Previous 
methods 


Turtle 
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Under a microscope, the filaments appear as a sequence of cylin- 
ders of various lengths, with a-type cells longer than b-type cells. The 
corresponding schematic image of filament development is shown in 
Figure 1.4. Note that due to the discrete nature of L-systems, the con- 
tinuous growth of cells between subdivisions is not captured by this 
model. 
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The geometric interpretation of strings applied to generate schematic 
images of Anabaena catenula is a very simple one. Letters of the 
L-system alphabet are represented graphically as shorter or longer rect- 
angles with rounded corners. The generated structures are one-dimen- 
sional chains of rectangles, reflecting the sequence of symbols in the 
corresponding strings. 

In order to model higher plants, a more sophisticated graphical in- 
terpretation of L-systems is needed. The first results in this direction 
were published in 1974 by Frijters and Lindenmayer [46], and Hogeweg 
and Hesper [64]. In both cases, L-systems were used primarily to de- 
termine the branching topology of the modeled plants. The geometric 
aspects, such as the lengths of line segments and the angle values, were 
added in a post-processing phase. The results of Hogeweg and Hesper 
were subsequently extended by Smith [136, 137], who demonstrated the 
potential of L-systems for realistic image synthesis. 

Szilard and Quinton [141] proposed a different approach to L-system 
interpretation in 1979. They concentrated on image representations 
with rigorously defined geometry, such as chain coding [43], and showed 
that strikingly simple DOL-systems could generate the intriguing, con- 
voluted curves known today as fractals [95]. These results were sub- 
sequently extended in several directions. Siromoney and Subrama- 
nian [135] specified L-systems which generate classic space-filling curves. 
Dekking investigated the limit properties of curves generated by L- 
systems [32] and concentrated on the problem of determining the fractal 
(Hausdorff) dimension of the limit set [31]. Prusinkiewicz focused on 
an interpretation based on a LOGO-style turtle [1] and presented more 
examples of fractals and plant-like structures modeled using L-systems 
[109, 111]. Further applications of L-systems with turtle interpretation 
include realistic modeling of herbaceous plants [117], description of ko- 
lam patterns (an art form from Southern India) [112, 115, 133, 134], 
synthesis of musical scores [110] and automatic generation of space- 
filling curves [116]. 

The basic idea of turtle interpretation is given below. A state of the 
turtle is defined as a triplet (x, y,a@), where the Cartesian coordinates 
(x, y) represent the turtle’s position, and the angle a, called the heading, 
is interpreted as the direction in which the turtle is facing. Given 
the step size d and the angle increment 6, the turtle can respond to 


1.3. Turtle interpretation of strings i 


a b 


Start 
FFF-FF-F-F+F+FF-F-FFF 


Figure 1.5: (a) Turtle interpretation of string symbols F’, +, —. (b) Inter- 
pretation of a string. The angle increment 6 is equal to 90°. Initially the 
turtle faces up. 


commands represented by the following symbols (Figure 1.5a): 


F Move forward a step of length d. The state of the turtle 
changes to (2’,y’,a@), where z’ = x +dcosa and y’ = 
y+dsina. A line segment between points (x,y) and 
(x’, y’) is drawn. 


f Move forward a step of length d without drawing a line. 


+ Turn left by angle 6. The next state of the turtle is 
(x,y,a+06). The positive orientation of angles is counter- 
clockwise. 


— Turn right by angle 6. The next state of the turtle is 
(x, Y,a— 0). 


Given a string v, the initial state of the turtle (xo, yo, ao) and fixed 
parameters d and 6, the turtle interpretation of v is the figure (set of 
lines) drawn by the turtle in response to the string v (Figure 1.5b). 
Specifically, this method can be applied to interpret strings which are 
generated by L-systems. For example, Figure 1.6 presents four approxi- 
mations of the quadratic Koch island taken from Mandelbrot’s book |95, 
page 51]. These figures were obtained by interpreting strings generated 
by the following L-system: 


io hahha 
pp: PSPSPLPL PPE PHP SE 


The images correspond to the strings obtained in derivations of length 
0 to 3. The angle increment 6 is equal to 90°. The step length d is 
decreased four times between subsequent images, making the distance 


Interpretation 
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a b 
d 


oe 
22 BP 


Figure 1.6: Generating a quadratic Koch island 


between the endpoints of the successor polygon equal to the length of 
the predecessor segment. 

The above example reveals a close relationship between Koch con- 
structions and L-systems. The initiator corresponds to the axiom and 
the generator corresponds to the production successor. The predeces- 
sor F' represents a single edge. L-systems specified in this way can be 
perceived as codings for Koch constructions. Figure 1.7 presents further 
examples of Koch curves generated using L-systems. A slight compli- 
cation occurs if the curve is not connected; a second production (with 
the predecessor f) is then required to keep components the proper dis- 
tance from each other (Figure 1.8). The ease of modifying L-systems 
makes them suitable for developing new Koch curves. For example, one 
can start from a particular L-system and observe the results of insert- 
ing, deleting or replacing some symbols. A variety of curves obtained 
this way are shown in Figure 1.9. 


Koch 
constructions 
us. L-systems 
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n=2, 6=90° n=4, 6=90° 
F-F-F-F -F 
F — F+FF-FF-F-F+F+F F — F+F-F-F+F 


F-F-F+F+FF+FF-F 


Figure 1.7: Examples of Koch curves generated using L-systems: (a) 
Quadratic Koch island [95, page 52], (b) A quadratic modification of the 
snowflake curve [95, page 139] 





n=2, 6=90° 

F+F+F+F 

F > F+f-FF+F+FF+F6+FF-f+FF-F-FF-FEFFF 
f — ffffff 


Figure 1.8: Combination of islands and lakes [95, page 121] 
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= 90° 
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n=3, 6 
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n=5, 6=90° 
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e€ 





F 


F-F 


-F-F 


F+F 


FF 


F 


F > F-FF--F 


Figure 1.9: A sequence of Koch curves obtained by successive modification 


of the production successor 
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a n=10, 5=90° b n=6, 5=60° 
Fy Fy 
Fi, —-F)+F,+ Fi —-F,+F,+F, 
F,—-F,-F, F,—F,-F,-Fy 


Figure 1.10: Examples of curves generated by edge-rewriting L-systems: (a) 
the dragon curve [48], (b) the Sierpiriski gasket [132] 


1.4 Synthesis of DOL-systems 


Random modification of productions gives little insight into the rela- 
tionship between L-systems and the figures they generate. However, 
we often wish to construct an L-system which captures a given struc- 
ture or sequence of structures representing a developmental process. 
This is called the inference problem in the theory of L-systems. Al- 
though some algorithms for solving it were reported in the literature 
[79, 88, 89], they are still too limited to be of practical value in the 
modeling of higher plants. Consequently, the methods introduced be- 
low are more intuitive in nature. They exploit two modes of operation 
for L-systems with turtle interpretation, called edge rewriting and node 
rewriting using terminology borrowed from graph grammars [56, 57, 87]. 
In the case of edge rewriting, productions substitute figures for poly- 
gon edges, while in node rewriting, productions operate on polygon 
vertices. Both approaches rely on capturing the recursive structure of 
figures and relating it to a tiling of a plane. Although the concepts are 
illustrated using abstract curves, they apply to branching structures 
found in plants as well. 


1.4.1 Edge rewriting 


Edge rewriting can be viewed as an extension of Koch constructions. 
For example, Figure 1.10a shows the dragon curve [21, 48, 95] and the 
L-system that generated it. Both the F; and F, symbols represent 
edges created by the turtle executing the “move forward” command. 
The productions substitute F, or F, edges by pairs of lines forming 
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Figure 1.11: Examples of FASS curves generated by edge-rewriting L- 
systems: (a) hexagonal Gosper curve [51], (b) quadratic Gosper curve [32] 
or E-curve [96] 


left or right turns. Many interesting curves can be obtained assuming 
two types of edges, “left” and “right.” Figures 1.10b and 1.11 present 
additional examples. 

The curves included in Figure 1.11 belong to the class of FASS 
curves (an acronym for space-filling, self-avoiding, simple and self- 
similar) [116], which can be thought of as finite, self-avoiding approxi- 
mations of curves that pass through all points of a square (space-filling 
curves [106]). McKenna [96] presented an algorithm for constructing 
FASS curves using edge replacement. It exploits the relationship be- 
tween such a curve and a recursive subdivision of a square into tiles. 
For example, Figure 1.12 shows the tiling that corresponds to the E- 
curve of Figure 1.11b. The polygon replacing an edge F; (Figure 1.12a) 
approximately fills the square on the left side of F; (b). Similarly, the 
polygon replacing an edge F;. (c) approximately fills the square on the 
right side of that edge (d). Consequently, in the next derivation step, 
each of the 25 tiles associated with the curves (b) or (d) will be covered 
by their reduced copies (Figure 1.11b). A recursive application of this 
argument indicates that the whole curve is approximately space-filling. 
It is also self-avoiding due to the following two properties: 
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FF, F,+F,+F,-F,-F)+F,+F,F)-F,-Fi FF + 
F,-F,-F)F,-F,+F F,+F,+Fy-F1-F,F,+ 

F,—-F]F,+F,+F,-F]-F)F,-Fi+F,F,+FitF,- 
FFF, +F, +F,F)-F)-F,+F,+F1-F)-F,F; 


Figure 1.12: Construction of the E-curve on the square grid. Left and right 
edges are distinguished by the direction of ticks. 


e the generating polygon is self-avoiding, and 


e no matter what the relative orientation of the polygons lying on 
two adjacent tiles, their union is a self-avoiding curve. 


The first property is obvious, while the second can be verified by con- 
sidering all possible relative positions of a pair of adjacent tiles. 

Using a computer program to search the space of generating poly- 
gons, McKenna found that the E-curve is the simplest FASS curve 
obtained by edge replacement in a square grid. Other curves require 
generators with more edges (Figure 1.13). The relationship between 
edge rewriting and tiling of the plane extends to branching structures, 
providing a method for constructing and analyzing L-systems which 
operate according to the edge-rewriting paradigm (see Section 1.10.3). 


1.4.2 Node rewriting 


The idea of node rewriting is to substitute new polygons for nodes of the 
predecessor curve. In order to make this possible, turtle interpretation 
is extended by symbols which represent arbitrary subfigures. As shown 
in Figure 1.14, each subfigure A from a set of subfigures A is represented 
by: 


e two contact points, called the entry point P4 and the exit point 
Qa, and 


e two direction vectors, called the entry vector p,4 and the exit vector 
qa- 

During turtle interpretation of a string v, asymbol A € A incorporates 

the corresponding subfigure into a picture. To this end, A is translated 
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edge replacement: (a) a SquaRecurve (grid size 7 x 7), (b) an E-tour (grid 


Figure 1.13: Examples of FASS curves generated on the square grid using 
size 9 x 9). Both curves are from [96]. 





Figure 1.14: Description of a subfigure A 
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Figure 1.15: Recursive construction of the Hilbert curve [63] in terms of 
node replacement 


and rotated in order to align its entry point P, and direction p’4 with 
the current position and orientation of the turtle. Having placed A, the 
turtle is assigned the resulting exit position Q4 and direction (4. 

For example, assuming that the contact points and directions of 
subfigures L,, and R,, are as in Figure 1.15, the figures L,,,,; and R,4+1 
are captured by the following formulas: 


tRyaF = 1g k ly =F iat 
=L, + A,P A, > Pi 


Ling 
n+1 


Suppose that curves Lo and Ro are given. One way of evaluating the 
string Ly, (or iy) for n > 0 is to generate successive strings recur- 
sively, in the order of decreasing value of index n. For example, the 
computation of Lz would proceed as follows: 


is 2 SR SA 
S (oP yh Rg Lg S CLRoP Soh lg= PR 
FU ERo = Lok lg FR) SF lof eho Flip 





Recursive 
formulas 


Pure curves 
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Thus, the generation of string L,, can be considered as a string-rewriting 
mechanism, where the symbols on the left side of the recursive formulas 
are substituted by corresponding strings on the right side. The substi- 
tution proceeds in a parallel way with F’,+ and — replacing themselves. 
Since all indices in any given string have the same value, they can be 
dropped, provided that a global count of derivation steps is kept. Con- 
sequently, string L, can be obtained in a derivation of length n using 
the following L-system: 


w: LD 
py: Lo+RF-LFL- FR+ 
po: R-o-LF+RFR+FL- 


In order to complete the curve definition, it is necessary to specify 
the subfigures represented by symbols L and R. In the special case of 
pure curves [116], these subfigures are reduced to single points. Thus, 
one can assume that symbols L and R are erased (replaced by the empty 
string) at the end of the derivation. Alternatively, they can be left in 
the string and ignored by the turtle during string interpretation. This 
second approach is consistent with previous definitions of turtle inter- 
pretation [109, 112]. A general discussion of the relationship between 
recurrent formulas and L-systems is presented in [61, 62]. 

Construction of the L-system generating the Hilbert curve can be 
extended to other FASS curves [116]. Consider an array of m xm square 
tiles, each including a smaller square, called a frame. The edges of the 
frame run at some distance from the tile’s edges. Each frame bounds an 
open self-avoiding polygon. The endpoints of this polygon coincide with 
the two contact vertices of the frame. Suppose that a single-stroke line 
running through all tiles can be constructed by connecting the contact 
vertices of neighboring frames using short horizontal or vertical line 
segments. A FASS curve can be constructed by the recursive repetition 
of this connecting pattern. To this end, the array of m x m connected 
tiles is considered a macrotile which contains an open polygon inscribed 
into a macroframe. An array of m x m macrotiles is formed, and the 
polygons inscribed into the macroframes are connected together. This 
construction is carried out recursively, with m x m macrotiles at level 
n yielding one macrotile at level n + 1. 

Tile arrangements suitable for the generation of FASS curves can 
be found algorithmically, by searching the space of all possible arrange- 
ments on a grid of a given size. Examples of curves synthesized this 
way are given in Figures 1.16 and 1.17. 

As in the case of edge rewriting, the relationship between node 
rewriting and tilings of the plane extends to branching structures. It 
offers a method for synthesizing L-systems that generate objects with a 
given recursive structure, and links methods for plant generation based 
on L-systems with those using iterated function systems [7] (see Chap- 
ter 8). 
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b 
a n=2, d=90° 
n=3 , d=90° -L 
-L L—LFLF+RFR+FLFL-FRF-LFL- 
L—LF+RFR+FL-F-LFLFL-FRFR+ FR+F+RF-LFL-FRFRFR+ 


R—-LFLF+RFRFR+F+RF-LFL-FR R—-LFLFLF+RFR+FL-F-LF+RFR+ 
FLF+RFRF-LFL-FRFR 


Figure 1.16: Sample FASS curves constructed using tiles with contact points 
positioned along a tile edge: (a) 3 x 3 tiles form a macrotile, (b) 4 x 4 tiles 
form a macrotile 


a b 


a n=2, 0=90° 
1 
L—LFRFL-F-RFLFR+F+LFRFL 
R—RFLFR+F+LFRFL-F-RFLFR 


b n=2, 6=45° 
L 
L—L+F+R-F-L+F+R-F-L-F-R+F+L-F-R-F-L+F+R-F-L-F-R-F- 
L+F+R+F+L+F+R-F-L+F+R+F+L-R-F+F+L+F+R-F-L+F+R-F-L 
R—R-F-L+F+R-F-L+F+R+F+L-F-R+F+L+F+R-F-L+F+R+F+L+F+ 
R-F-L-F-R-F-L+F+R-F-L-F-R+F+L-F-R-F-L+F+R-F-L+F+R 


Figure 1.17: Sample FASS curves constructed using tiles with contact points 
positioned diagonally: (a) 3 x 3 tiles form a macrotile (Peano curve [106]), 
(b) 5 x 5 tiles form a macrotile 
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1.4.3 Relationship between edge and 
node rewriting 


The classes of curves that can be generated using the edge-rewriting and 
node-rewriting techniques are not disjoint. For example, reconsider the 
L-system that generates the dragon curve using edge replacement: 


w: FF 
pis fia ee 
P2: Lp Shy Ee 


Assume temporarily that a production predecessor can contain more 
than one letter; thus an entire subword can be replaced by the successor 
of a single production (a formalization of this concept is termed pseudo- 
L-systems {109]). The dragon-generating L-system can be rewritten as: 


w: Fl 
pp: Flo Fl+rF+ 
po: rE > —-Fl—rF 


where the symbols / and r are not interpreted by the turtle. Production 
p1 replaces the letter / by the string 1+ rF— while the leading letter F’ 
is left intact. In a similar way, production pz replaces the letter r by 
the string —F'l—r and leaves the trailing F intact. Thus, the L-system 
can be transformed into node-rewriting form as follows: 


w: Fl 
py: lol+rF+ 
055 CS SSF 


In practice, the choice between edge rewriting and node rewriting 
is often a matter of convenience. Neither approach offers an auto- 
matic, general method for constructing L-systems that capture given 
structures. However, the distinction between edge and node rewriting 
makes it easier to understand the intricacies of L-system operation, and 
in this sense assists in the modeling task. Specifically, the problem of 
filling a region by a self-avoiding curve is biologically relevant, since 
some plant structures, such as leaves, may tend to fill a plane without 
overlapping [38, 66, 67, 94]. 


1.5 Modeling in three dimensions 


Turtle interpretation of L-systems can be extended to three dimensions 
following the ideas of Abelson and diSessa [1]. The key concept is to 
represent the current orientation of the turtle in space by three vectors 
H,L,U, indicating the turtle’s heading, the direction to the left, and the 
direction up. These vectors have unit length, are perpendicular to each 


1.5. Modeling in three dimensions 19 





Figure 1.18: Controlling the turtle in three dimensions 


other, and satisfy the equation # x Fr = U. Rotations of the turtle are 
then expressed by the equation 


(A i o)=[a z a]R, 
where R is a 3 x 3 rotation matrix [40]. Specifically, rotations by angle 
a about vectors U,£ and dH are represented by the matrices: 
cosa sina 0 


Ru(a) = | —sina cosa 0 
0 0 1 


cosa 0 —sina 


Rr (a) = | 0 1 0 | 


sina 0 cosa 


I.-Y 0 
Ru(a) = | 0 cosa —sina 

0 sina cosa 
The following symbols control turtle orientation in space (Figure 1.18): 
-- Turn left by angle 6, using rotation matrix Ry(0). 
- Turn right by angle 6, using rotation matrix Ry(—4). 
Pitch down by angle 6, using rotation matrix Ry(0). 
Pitch up by angle 6, using rotation matrix Ry(—9). 


Roll left by angle 6, using rotation matrix Ry(0). 


Se. wee” re RE 


Roll right by angle 6, using rotation matrix Ry(—6). 


| Turn around, using rotation matrix Ry(180°). 
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B-F+CFC+F-D&F AD-F+&&CFC+F+B// 
A&F\CFBAFADAA-F-DA|FAB|FCAFAA// 
|DA|FAB-F+CAFAA&&FA&FAC+F+BAFAD// 
|CFB-F+B|FA&F\A&&FB-F+B|FC// 


Figure 1.19: A three-dimensional extension of the Hilbert curve [139]. Col- 
ors represent three-dimensional “frames” associated with symbols A (red), B 
(blue), C (green) and D (yellow). 
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As an example of a three-dimensional object created using an L- 
system, consider the extension of the Hilbert curve shown in Figure 1.19. 
The L-system was constructed with the node-replacement technique 
discussed in the previous section, using cubes and “macrocubes” in- 
stead of tiles and macrotiles. 


1.6 Branching structures 


According to the rules presented so far, the turtle interprets a character 
string as a sequence of line segments. Depending on the segment lengths 
and the angles between them, the resulting line is self-intersecting or 
not, can be more or less convoluted, and may have some segments 
drawn many times and others made invisible, but it always remains 
just a single line. However, the plant kingdom is dominated by branch- 
ing structures; thus a mathematical description of tree-like shapes and 
methods for generating them are needed for modeling purposes. An 
axial tree [89, 117] complements the graph-theoretic notion of a rooted 
tree [108] with the botanically motivated notion of branch axis. 


1.6.1 Axial trees 


A rooted tree has edges that are labeled and directed. The edge se- 
quences form paths from a distinguished node, called the root or base, 
to the terminal nodes. In the biological context, these edges are re- 
ferred to as branch segments. A segment followed by at least one more 
segment in some path is called an internode. A terminal segment (with 
no succeeding edges) is called an apez. 

An avzial tree is a special type of rooted tree (Figure 1.20). At each 
of its nodes, at most one outgoing straight segment is distinguished. 
All remaining edges are called lateral or side segments. A sequence of 
segments is called an azis if: 


e the first segment in the sequence originates at the root of the tree 
or as a lateral segment at some node, 


e each subsequent segment is a straight segment, and 


e the last segment is not followed by any straight segment in the 
tree. 


Together with all its descendants, an axis constitutes a branch. A 
branch is itself an axial (sub)tree. 

Axes and branches are ordered. The axis originating at the root of 
the entire plant has order zero. An axis originating as a lateral segment 
of an n-order parent axis has order n+1. The order of a branch is equal 
to the order of its lowest-order or main axis. 
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Figure 1.20: An axial tree 





Figure 1.21: 


Sample tree generated using a method based on Horton— 


Strahler analysis of branching patterns 
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Figure 1.22: A tree production p and its application to the edge S in a tree 
Ty 


Axial trees are purely topological objects. The geometric connotation 
of such terms as straight segment, lateral segment and axis should be 
viewed at this point as an intuitive link between the graph-theoretic 
formalism and real plant structures. 

The proposed scheme for ordering branches in axial trees was in- 
troduced originally by Gravelius [53]. MacDonald [94, pages 110-121] 
surveys this and other methods applicable to biological and geograph- 
ical data such as stream networks. Of special interest are methods 
proposed by Horton [70, 71] and Strahler, which served as a basis for 
synthesizing botanical trees [37, 152] (Figure 1.21). 


1.6.2 Tree OL-systems 


In order to model development of branching structures, a rewriting 
mechanism can be used that operates directly on axial trees. A rewrit- 
ing rule, or tree production, replaces a predecessor edge by a successor 
axial tree in such a way that the starting node of the predecessor is 
identified with the successor’s base and the ending node is identified 
with the successor’s top (Figure 1.22). 

A tree OL-system G is specified by three components: a set of edge 
labels V, an initial tree w with labels from V, and a set of tree produc- 
tions P. Given the L-system G, an axial tree 7) is directly derived from 
a tree T,, noted 7, => To, if T> is obtained from 7 by simultaneously 
replacing each edge in 7) by its successor according to the production 
set P. A tree T is generated by G in a derivation of length n if there 
exists a sequence of trees To, 7;,...,7;, such that To = w, 7, = T and 
Tp =>T,=>...=> Tj. 


Stack 
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Zz, 


A F[+F] [-F[-F] F] F[+F] [-F] 


Figure 1.23: Bracketed string representation of an axial tree 


1.6.3 Bracketed OL-systems 


The definition of tree L-systems does not specify the data structure for 
representing axial trees. One possibility is to use a list representation 
with a tree topology. Alternatively, axial trees can be represented using 
strings with brackets [82]. A similar distinction can be observed in Koch 
constructions, which can be implemented either by rewriting edges and 
polygons or their string representations. An extension of turtle in- 
terpretation to strings with brackets and the operation of bracketed 
L-systems [109, 111] are described below. 

Two new symbols are introduced to delimit a branch. They are 
interpreted by the turtle as follows: 


[ Push the current state of the turtle onto a pushdown 
stack. The information saved on the stack contains the 
turtle’s position and orientation, and possibly other at- 
tributes such as the color and width of lines being drawn. 


| Pop a state from the stack and make it the current state 
of the turtle. No line is drawn, although in general the 
position of the turtle changes. 


An example of an axial tree and its string representation are shown 
in Figure 1.23. 

Derivations in bracketed OL-systems proceed as in OL-systems with- 
out brackets. The brackets replace themselves. Examples of two- 
dimensional branching structures generated by bracketed OL-systems 
are shown in Figure 1.24. 

Figure 1.25 is an example of a three-dimensional bush-like structure 
generated by a bracketed L-system. Production p; creates three new 
branches from an apex of the old branch. A branch consists of an 
edge F’ forming the initial internode, a leaf L and an apex A (which 
will subsequently create three new branches). Productions pz and ps3 
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a b c 
n=5 ,0=25.7° n=5 ,d=20° ere 

F F - 
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[+F-F-F] 








d e f 

n=7 ,0=20° n=7 ,0=25.7° n=5 ,0=22.5° 

X X X 

X —F [+X] F [-X]+X X—F [+X] [-X] FX X—F- [[X]+X]+F [+FX] -X 
F FF F—FF F—FF 


Figure 1.24: Examples of plant-like structures generated by bracketed OL- 
systems. L-systems (a), (b) and (c) are edge-rewriting, while (d), (e) and 
(f) are node-rewriting. 


Plant 
with flowers 


26 Chapter 1. Graphical modeling using L-systems 





n=7, 0=22.5° 

Ww A 

pi: A — [&FL!IA)/////? [(&FL!A]///////? (&FL!A) 
po: FS ///// F 

pa: SFL 

pa Lem PPOPAAL HE tee ee} | 


Figure 1.25: A three-dimensional bush-like structure 


specify internode growth. In subsequent derivation steps the internode 
gets longer and acquires new leaves. This violates a biological rule 
of subapical growth (discussed in detail in Chapter 3), but produces 
an acceptable visual effect in a still picture. Production p, specifies 
the leaf as a filled polygon with six edges. Its boundary is formed 
from the edges f enclosed between the braces { and } (see Chapter 5 
for further discussion). The symbols ! and ' are used to decrement 
the diameter of segments and increment the current index to the color 
table, respectively. 

Another example of a three-dimensional plant is shown in Fig- 
ure 1.26. The L-system can be described and analyzed in a way similar 
to the previous one. 
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n=5, 6=18° 
w : plant 
pi: plant — internode + [ plant + flower] — — // 
[ — — leaf | internode | + + leaf ] — 
[ plant flower | + + plant flower 
p2: internode — F seg [// & & leaf | [// A A leaf | F seg 
p3: seg — seg F seg 
pa: leaf [ { +f-ff-f+ | +f-ff-f } 
ps: flower > [& & & pedicel ‘ / wedge //// wedge //// 
wedge //// wedge //// wedge | 
pe: pedicel — FF 
pr: wedge>['AF][{&&& & —f4+f | —f4+f } | 


Figure 1.26: A plant generated by an L-system 
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1.7 Stochastic L-systems 


All plants generated by the same deterministic L-system are identical. 
An attempt to combine them in the same picture would produce a 
striking, artificial regularity. In order to prevent this effect, it is nec- 
essary to introduce specimen-to-specimen variations that will preserve 
the general aspects of a plant but will modify its details. 

Variation can be achieved by randomizing the turtle interpretation, 
the L-system, or both. Randomization of the interpretation alone has 
a limited effect. While the geometric aspects of a plant — such as 
the stem lengths and branching angles — are modified, the underly- 
ing topology remains unchanged. In contrast, stochastic application 
of productions may affect both the topology and the geometry of the 
plant. The following definition is similar to that of Yokomori [162] and 
Eichhorst and Savitch [35]. 

A stochastic OL-system is an ordered quadruplet G, = (V,w,P,7). 
The alphabet V, the axiom w and the set of productions P are defined 
as in an OL-system (page 4). Function 7 : P — (0,1), called the 
probability distribution, maps the set of productions into the set of 
production probabilities. It is assumed that for any letter a € V, the 
sum of probabilities of all productions with the predecessor a is equal 
to 1. 

We will call the derivation pp => v a stochastic derivation in G, if for 
each occurrence of the letter a in the word yz the probability of applying 
production p with the predecessor a is equal to m(p). Thus, different 
productions with the same predecessor can be applied to various occur- 
rences of the same letter in one derivation step. 

A simple example of a stochastic L-system is given below. 





wi 

pp: F 3S FU+F\F[-FI|F 
po. FS FI+F\F 

p3. FS F[-F\F 


The production probabilities are listed above the derivation symbol 
—. Each production can be selected with approximately the same 
probability of 1/3. Examples of branching structures generated by this 
L-system with derivations of length 5 are shown in Figure 1.27. Note 
that these structures look like different specimens of the same (albeit 
fictitious) plant species. 

A more complex example is shown in Figure 1.28. The field consists 
of four rows and four columns of plants. All plants are generated by a 
stochastic modification of the L-system used to generate Figure 1.26. 
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Figure 1.28: Flower field 
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The essence of this modification is to replace the original production 
p3 by the following three productions: 


ph: seg = seg [ // & & leaf] [// AA leaf | F seg 


py: seg 23 seg F seg 


Mm, 34 
py: seg > seg 


Thus, in any step of the derivation, the stem segment seg may either 
grow and produce new leaves (production p3), grow without producing 
new leaves (production p3), or not grow at all (production p4’). All three 
events occur with approximately the same probability. The resulting 
field appears to consist of various specimens of the same plant species. 
If the same L-system was used again (with different seed values for the 
random number generator), a variation of this image would be obtained. 


1.8 Context-sensitive L-systems 


Productions in OL-systems are context-free; i.e. applicable regardless 
of the context in which the predecessor appears. However, production 
application may also depend on the predecessor’s context. This effect is 
useful in simulating interactions between plant parts, due for example to 
the flow of nutrients or hormones. Various context-sensitive extensions 
of L-systems have been proposed and studied thoroughly in the past 
(62, 90, 128]. 2L-systems use productions of the form a; < a >a, — x, 
where the letter a (called the strict predecessor) can produce word y if 
and only if a is preceded by letter a; and followed by a,. Thus, letters 
a, and a, form the left and the right context of a in this production. 
Productions in 1L-systems have one-sided context only; consequently, 
they are either of the form a; < a — y ora > a, — y. OL-systems, 
1L-systems and 2L-systems belong to a wider class of [L-systems, also 
called (k,l)-systems. In a (k,l)-system, the left context is a word of 
length k and the right context is a word of length J. 

In order to keep specifications of L-systems short, the usual notion 
of IL-systems has been modified here by allowing productions with 
different context lengths to coexist within a single system. Further- 
more, context-sensitive productions are assumed to have precedence 
over context-free productions with the same strict predecessor. Conse- 
quently, if a context-free and a context-sensitive production both apply 
to a given letter, the context-sensitive one should be selected. If no pro- 
duction applies, this letter is replaced by itself as previously assumed 
for OL-systems. 
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Figure 1.29: The predecessor of a context-sensitive tree production (a) 
matches edge S in a tree T (b) 


The following sample 1L-system makes use of context to simulate signal Signal 
propagation throughout a string of symbols: propagation 


w : baaaaaaaa 
pi: b<a —b 
Do: b a 


The first few words generated by this L-system are given below: 


baaaaaaaa 
abaaaaaaa 
aabaaaaaa 
aaabaaaaa 
aaaabaaaa 


The letter b moves from the left side to the right side of the string. 

A context-sensitive extension of tree L-systems requires neighbor Context in tree 
edges of the replaced edge to be tested for context matching. A prede- L-systems 
cessor of a context-sensitive production p consists of three components: 

a path | forming the left context, an edge S called the strict predecessor, 
and an axial tree r constituting the right context (Figure 1.29). The 
asymmetry between the left context and the right context reflects the 
fact that there is only one path from the root of a tree to a given edge, 
while there can be many paths from this edge to various terminal nodes. 
Production p matches a given occurrence of the edge S in a tree T if 1 
is a path in T terminating at the starting node of S, and r is a subtree 


Context in 
bracketed 
L-systems 
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of T originating at the ending node of S. The production can then be 
applied by replacing S with the axial tree specified as the production 
successor. 

The introduction of context to bracketed L-systems is more difficult 
than in L-systems without brackets, because the bracketed string repre- 
sentation of axial trees does not preserve segment neighborhood. Conse- 
quently, the context matching procedure may need to skip over symbols 
representing branches or branch portions. For example, Figure 1.29 in- 
dicates that a production with the predecessor BC < S > G|H]M can 
be applied to symbol S in the string 


ABC|DE|[SG[HI[.JK|L|M NO}, 


which involves skipping over symbols |D EF] in the search for left context, 
and I|JK|L in the search for right context. 

Within the formalism of bracketed L-systems, the left context can 
be used to simulate control signals that propagate acropetally, i.e., from 
the root or basal leaves towards the apices of the modeled plant, while 
the right context represents signals that propagate baszpetally, t.e., from 
the apices towards the root. For example, the following 1L-system 
simulates propagation of an acropetal signal in a branching structure 
that does not grow: 


#ignore : +— 


wi: Fil+F,|Fo|—F.]Po[+Fal Pa 
Pie. Peay Fy 


Symbol Ff; represents a segment already reached by the signal, while 
F, represents a segment that has not yet been reached. The #ignore 
statement indicates that the geometric symbols + and — should be 
considered as non-existent while context matching. Images representing 
consecutive stages of signal propagation (corresponding to consecutive 
words generated by the L-system under consideration) are shown in 
Figure 1.30a. 

The propagation of a basipetal signal can be simulated in a similar 
way (Figure 1.30b): 


#ignore : +— 


w: BGP IE RIFLE |F 
Pi: F, > F, _- F, 
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Figure 1.30: Signal propagation in a branching structure: (a) acropetal, (b) 
basipetal 


The operation of context-sensitive L-systems is examined further using 
examples obtained by Hogeweg and Hesper [64]. In 1974, they pub- 
lished the results of an exhaustive study of 3,584 patterns generated 
by a class of bracketed 2L-systems defined over the alphabet {0,1}. 
Some of these patterns had plant-like shapes. Subsequently, Smith 
significantly improved the quality of the generated images using state- 
of-the-art computer imagery techniques [136, 137]. Sample structures 
generated by L-systems similar to those proposed by Hogeweg and Hes- 
per are shown in Figure 1.31. The differences are related to the geo- 
metric interpretation of the resulting strings. According to the original 
interpretation, consecutive branches are issued alternately to the left 
and right, whereas turtle interpretation requires explicit specification 
of branching angles within the L-system. 
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Figure 1.31: Examples of branching structures generated using L-systems 


based on the results of Hogeweg and Hesper [64] 
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Figure 1.31 (continued): L-systems of Hogeweg and Hesper 
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1.9 Growth functions 


During the synthesis of a plant model it is often convenient to dis- 
tinguish productions that specify the branching pattern from those 
that describe elongation of plant segments. This separation can be 
observed in some of the L-systems considered so far. For example, in 
L-systems (d), (e) and (f) from Figure 1.24 the first productions cap- 
ture the branching patterns, while the remaining productions, equal in 
all cases to F — F'F, describe elongation of segments represented by 
sequences of symbols F’. The number of letters F' in a string vy, orig- 
inating from a single letter F' is doubled in each derivation step, thus 
the elongation is exponential, with length(y,,) = 2”. 

A function that describes the number of symbols in a word in terms 
of its derivation length is called a growth function. The theory of L- 
systems contains an extensive body of results on growth functions [62, 
127]. The central observation is that the growth functions of DOL- 
systems are independent of the letter ordering in the productions and 
derived words. Consequently, the relation between the number of letter 
occurrences in a pair of words yp and v, such that ~ => v, can be 
conveniently expressed using matrix notation. 

Let G = (V,w,P) be a DOL-system and assume that letters of 
the alphabet V have been ordered, V = {a1, d2,...,@m}. Construct a 
square matrix Qmxm, Where entry qj; is equal to the number of occur- 
rences of letter a; in the successor of the production with predecessor 
a;. Let a* denote the number of occurrences of letter a; in the word 
x generated by G in a derivation of length k. The definition of direct 
derivation in a DOL-system implies that 


dir G12 tt) Gm 
G21 G22 *** Gam 

[ af ak ..- ak, = fiaee aktl 1. ght 
Ami Um2 °*** Umm 


This matrix notation is useful in the analysis of growth functions. For 
example, consider the following L-system: 


to WG 
pi: a—ab (1.2) 
po: boa 


The relationship between the number of occurrences of letters a and b 
in two consecutively derived words is 


[ a ly ober pet! 
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or 
a®t! — gk + b* = a* + a 
tor k = 1,2,3,..... From the axiom it follows that a® = 1 and 
a = pf = 0. Thus, the number of letters a in the strings gener- 
ated by the L-system specified i in equation (1.2) grows according to the 
Fibonacci series: 1,1,2,3,5,8,... . This growth function was imple- 
mented by productions py and Da ‘in the L-system generating the bush 
in Figure 1.25 (page 26) to describe the elongation of its internodes. 
Polynomial growth functions of arbitrary degree can be obtained 
using L-systems of the following form: 





Ww ao 

Pi: Qo aoa, 
P2: G1 —~ a1a2 
P3: G2 — a2az 
Pa: G3 — A304 


The matrix Q is given below: 


1 

0 

asl 0 
- 0 


Eee OO 
| 


SGOrHE 
- ORFS 


Thus, for any i,k > 1, the number a? of occurrences of symbol a; in 
the string generated in a derivation of length k satisfies the equality 


ay + ai ny ait. 
Taking into consideration the axiom, the distribution of letters a; as 
a function of the derivation length is captured by the following table 


(only non-zero terms are shown): 


k k k k k k 
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|) ha? al 
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This table represents the Pascal triangle, thus for any k > i > 1 its 
terms satisfy the following equality: 


al oo 
= i Tedasah 

Consequently, the number of occurrences of letter a; as a function of 
the derivation length k is expressed by a polynomial of degree i. By 
identifying letter a; with the turtle symbol F’, it is possible to model in- 
ternode elongation expressed by polynomials of arbitrary degree 7 > 0. 
This observation was generalized by Szilard [140], who developed an al- 
gorithm for constructing a DOL-system with growth functions specified 
by any positive, nondecreasing polynomials with integer coefficients |62, 
page 276]. 

The examples of growth functions considered so far include expo- 
nential and polynomial functions. Rozenberg and Salomaa [127, pages 
30-38] show that, in general, the growth function fe(n) of any DOL- 
system G = (V,w, P) is a combination of polynomial and exponential 
functions: 


jot) = 3 P,(n)p; for n> no, (1.3) 
i=l 


where P;(n) denotes a polynomial with integer coefficients, p; is a non- 
negative integer, and no is the total number of letters in the alphabet 
V. Unfortunately, many growth processes observed in nature cannot be 
described by equation (1.3). Two approaches are then possible within 
the framework of the theory of L-systems. 

The first is to extend the size no of the alphabet V, so that the 
growth process of interest will be captured by the initial derivation 
steps, W = Lo => fly > +++ => Ln, before equation (1.3) starts to apply. 
For example, the L-system 


WwW : ao 
Pit A > i410 for 1<k (1.4) 
Pk+j : b; — Ona for j <l 


over the alphabet V = {ao, a1,..., @x}U {bo, bi, ..., bs } U{F'} can be used 
to approximate a sigmoidal elongation of a segment represented by a 
sequence of symbols F’ (Figure 1.32). The term sigmoidal refers to a 
function with a plot in the shape of the letter S. Such functions are 
commonly found in biological processes [143], with the initial part of 
the curve representing the growth of a young organism, and the latter 
part corresponding to the organism close to its final size. 

The second approach to the synthesis of growth functions out- 
side the class captured by equation (1.3) is to use context-sensitive 
L-systems. For example, the following 2L-system has a growth func- 
tion given by fe(n) = |./n| +4, where || is the floor function. 
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Figure 1.32: A sigmoidal growth function implemented using the L-system 
in equation (1.4), for k = 1 = 20 


w: XF,FLX 

mn: b< ho >h oO FF, 

po: Fy, <Fy>X 7 Faky 

p3: Fy<Fy>Fa - Fa (13) 
pa: X< > o KB 

Ps: Ee a Ee 

Pe: Fy = FS 


The operation of this L-system is illustrated in Figure 1.33. Produc- 
tions p; and p3, together with ps and pg, propagate symbols F), and Fy 
up and down the string of symbols jz. Productions pz and py change the 
propagation direction, after symbol X marking a string end has been 
reached by Fi, or Fy, respectively. In addition, pp extends the string 
with a symbol F,,. Thus, the number of derivation steps increases by 
two between consecutive applications of production po. As a result, 
string extension occurs at derivation steps n expressed by the square 
of the string length, which yields the growth function |./n] + 4. 

In practice it is often difficult, if not impossible, to find L-systems 
with the required growth functions. Vitdnyi [153] illustrates this by 
referring to sigmoidal curves: 


If we want to obtain sigmoidal growth curves with the 
original L-systems then not even the introduction of cell 
interaction can help us out. In the first place, we end up 
constructing quite unlikely flows of messages through the 
organism, which are more suitable to electronic computers, 
and in fact give the organism universal computing power. 
Secondly, and this is more fundamental, we can not obtain 
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Figure 1.33: Square-root growth implemented using the L-system specified 
in equation (1.5) 


growth which, always increasing the size of the organism, 
tends towards stability in the limit. The slowest increasing 
growth we can obtain by allowing cell interaction is loga- 
rithmic and thus can not at all account for the asymptotic 
behavior of sigmoidal growth functions. 


In the next section we present an extension of L-systems that makes 
it possible to avoid this problem by allowing for explicit inclusion of 
growth functions into L-system specifications. 


1.10 Parametric L-systems 


Although L-systems with turtle interpretation make it possible to gen- 
erate a variety of interesting objects, from abstract fractals to plant-like 
branching structures, their modeling power is quite limited. A major 
problem can be traced to the reduction of all lines to integer multi- 
ples of the unit segment. As a result, even such a simple figure as an 
isosceles right-angled triangle cannot be traced exactly, since the ratio 
of its hypotenuse length to the length of a side is expressed by the irra- 
tional number V2. Rational approximation of line length provides only 
a limited solution, because the unit step must be the smallest common 
denominator of all line lengths in the modeled structure. Consequently, 
the representation of a simple plant module, such as an internode, may 
require a large number of symbols. The same argument applies to an- 
gles. Problems become even more pronounced while simulating changes 
to the modeled structure over time, since some growth functions can- 
not be expressed conveniently using L-systems. Generally, it is difficult 
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to capture continuous phenomena, since the obvious technique of dis- 
cretizing continuous values may require a large number of quantization 
levels, yielding L-systems with hundreds of symbols and productions. 
Consequently, model specification becomes difficult, and the mathe- 
matical beauty of L-systems is lost. 

In order to solve similar problems, Lindenmayer proposed that nu- 
merical parameters be associated with L-system symbols [83]. He illus- 
trated this idea by referring to the continuous development of branching 
structures and diffusion of chemical compounds in a nonbranching fil- 
ament of Anabaena catenula. Both problems were revisited in later 
papers [25, 77]. A definition of parametric L-systems was formulated 
by Prusinkiewicz and Hanan [113] and is presented below. 


1.10.1 Parametric OL-systems 


Parametric L-systems operate on parametric words, which are strings 
of modules consisting of letters with associated parameters. The let- 
ters belong to an alphabet V, and the parameters belong to the set 
of real numbers ¥. A module with letter A € V and parameters 
1,42, -..,4n € Ris denoted by A(aj, az, ...,a,). Every module belongs 
to the set M = V x &*, where * is the set of all finite sequences of pa- 
rameters. The set of all strings of modules and the set of all nonempty 
strings are denoted by M* = (V x R*)* and M+ = (V x #*)t, respec- 
tively. 

The real-valued actual parameters appearing in the words corre- 
spond with formal parameters used in the specification of L-system 
productions. If © is a set of formal parameters, then C'() denotes a 
logical expression with parameters from /, and E(%) is an arithmetic 
expression with parameters from the same set. Both types of expres- 
sions consist of formal parameters and numeric constants, combined 
using the arithmetic operators +, —, *, /; the exponentiation operator 
A, the relational operators <, >, =; the logical operators !, &, | (not, 
and, or); and parentheses (). Standard rules for constructing syntac- 
tically correct expressions and for operator precedence are observed. 
Relational and logical expressions evaluate to zero for false and one for 
true. A logical statement specified as the empty string is assumed to 
have value one. The sets of all correctly constructed logical and arith- 
metic expressions with parameters from © are noted C(X) and €(%). 

A parametric OL-system is defined as an ordered quadruplet G = 
(V, %,w,P), where 


e V is the alphabet of the system, 
e » is the set of formal parameters, 
ew €(V x #*)* is a nonempty parametric word called the axiom, 


e PC(V x d*) x C(X) x (V x E(4))* is a finite set of productions. 
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The symbols : and — are used to separate the three components of a 
production: the predecessor, the condition and the successor. For exam- 
ple, a production with predecessor A(t), condition t > 5 and successor 
B(t + 1)CD(t A0.5,t — 2) is written as 


A(t) :t>5— B(t+1)CD(t A0.5,t — 2). (1.6) 


A production matches a module in a parametric word if the following 
conditions are met: 


e the letter in the module and the letter in the production prede- 
cessor are the same, 


e the number of actual parameters in the module is equal to the 
number of formal parameters in the production predecessor, and 


e the condition evaluates to true if the actual parameter values are 
substituted for the formal parameters in the production. 


A matching production can be applied to the module, creating a string 
of modules specified by the production successor. The actual parame- 
ter values are substituted for the formal parameters according to their 
position. For example, production (1.6) above matches a module A(9), 
since the letter A in the module is the same as in the production pre- 
decessor, there is one actual parameter in the module A(9) and one 
formal parameter in the predecessor A(t), and the logical expression 
t > 5 is true for t = 9. The result of the application of this production 
is a parametric word B(10)C'D(3, 7). 

If a module a produces a parametric word y as the result of a 
production application in an L-system G, we write a> y. Given a 
parametric word ~/ = @,@2...dm, we say that the word v = y1X2..-Xm 
is directly derived from (or generated by) js and write w => v if and 
only if a; > x; for all 2 = 1,2,...,.m. A parametric word v is generated 
by G in a derivation of length n if there exists a sequence of words 
Lo, #41) -++) ln Such that po = WwW, by = Vv and lo => py =... => Mn. 

An example of a parametric L-system is given below. 


w : B(2)A(4,4) 

pi: A(ayy) :y<=3 — Alxx2, x+y) 

po: A(z,y) :y>3 — B(x)A(z/y,0) (1.7) 
p3: Biz) :a@<1l - C 

pa: Birt) :x2>=1 > Bi(x-1) 


As in the case of non-parametric L-systems, it is assumed that a module 
replaces itself if no matching production is found in the set P. The 
words obtained in the first few derivation steps are shown in Figure 1.34. 
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Figure 1.34: The initial sequence of strings generated by the parametric 
L-system specified in equation (1.7) 


1.10.2 Parametric 2L-systems 


Productions in parametric OL-systems are context-free, 7.e., applicable 
regardless of the context in which the predecessor appears. A context- 
sensitive extension is necessary to model information exchange between 
neighboring modules. In the parametric case, each component of the 
production predecessor (the left context, the strict predecessor and the 
right context) is a parametric word with letters from the alphabet V 
and formal parameters from the set ©. Any formal parameters may 
appear in the condition and the production successor. 
A sample context-sensitive production is given below: 


A(z) < Bly) > C(z):e¢+y+2> 10> E((x+y)/2)F((y + z)/2) 
It can be applied to the module B(5) that appears in a parametric word 
-- A(4)B(5)C(6) --- (1.8) 


since the sequence of letters A, B, C in the production and in parametric 
word (1.8) are the same, the numbers of formal parameters and actual 
parameters coincide, and the condition 4+ 5+6 > 10 is true. Asa 
result of the production application, the module B(5) will be replaced 
by a pair of modules E(4.5)F(5.5). Naturally, the modules A(4) and 
C(6) will be replaced by other productions in the same derivation step. 

Parametric 2L-systems provide a convenient tool for expressing de- 
velopmental models that involve diffusion of substances throughout an 
organism. A good example is provided by an extended model of the 
pattern of cells observed in Anabaena catenula and other blue-green 
bacteria [99]. This model was proposed by de Koster and Linden- 
mayer [25]. 
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#define CH 900 /* high concentration */ 

#define CT 0.4 /* concentration threshold */ 
#define ST 3.9 /* segment size threshold */ 
#include H /* heterocyst shape specification */ 
#ignore f ~ H 


w : -(90)F(0,0,CH)F(4,1,CH)F(0O,0,CH) 

pi : F(s,t,c) : t=1 & s>=6 — 
F(s/3*2,2,c)fC1)FGs/3,.1,¢) 

po : F(s,t,c) : t=2 & s>=6 — 
F(s/3,2,c)f£(1)F(s/3*2,1,c) 

p3 : F(h,i,k) < F(s,t,c) > Flo,p,r) : s>STlc>CT -— 
F(st+.1,t,c+0.25* (k+r-3*c) ) 


pa : F(h,i,k) < F(s,t,c) > Flo,p,r) : !(s>STlc>CT) — 
F(0,0,CH) ~ H(1) 


ps : Hs) : s<3 — H(s*1.1) 


L-system 1.1: Anabaena catenula 


Generally, the bacteria under consideration form a nonbranching fila- 
ment consisting of two classes of cells: vegetative cells and heterocysts. 
Usually, the vegetative cells divide and produce two daughter vegeta- 
tive cells. This mechanism is captured by the L-system specified in 
equation (1.1) and Figure 1.4 (page 5). However, in some cases the 
vegetative cells differentiate into heterocysts. Their distribution forms 
a well-defined pattern, characterized by a relatively constant number 
of vegetative cells separating consecutive heterocysts. How does the 
organism maintain constant spacing of heterocysts while growing? The 
model explains this phenomenon using a biologically well-motivated 
hypothesis that heterocyst distribution is regulated by nitrogen com- 
pounds produced by the heterocysts, transported from cell to cell across 
the filament, and decayed in the vegetative cells. If the compound’s 
concentration in a young vegetative cell falls below a specific level, this 
cell differentiates into a heterocyst (L-system 1.1). 

The #define statements assign values to numerical constants used 
in the L-system. The #include statement specifies the shape of a het- 
erocyst (a disk) by referring to a library of predefined shapes (see Sec- 
tion 5.1). Cells are represented by modules F'(s,t,c), where s stands 
for cell length, ¢ is cell type (0 - heterocyst, 1 and 2 - vegetative types’), 





'The model of Anabaena introduced in Section 1.2 distinguished between four 
types of cells: a,, by, aj and b;. Cells b do not divide and can be considered as young 
forms of the corresponding cells a. Thus, the vegetative type 1 considered here 
embraces cells a, and b;, while type 2 embraces a; and b;. The formal relationship 
between the four-cell and two-cell models is further discussed in Chapter 6. 
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Figure 1.35: Development of Anabaena catenula with heterocysts, simulated 
using parametric L-system 1.1 


and c represents the concentration of nitrogen compounds. Productions 
p, and pz describe division of the vegetative cells. They each create two 
daughter cells of unequal length. The difference between cells of type 
1 and 2 lies in the ordering of the long and short daughter cells. Pro- 
duction p3 captures the processes of transportation and decay of the 
nitrogen compounds. Their concentration c is related to the concen- 
tration in the neighboring cells k and r, and changes in each derivation 
step according to the formula 


Cc =c+0.25(k+r—3 0). 


Production pq describes differentiation of a vegetative cell into a hete- 
rocyst. The condition specifies that this process occurs when the cell 
length does not exceed the threshold value ST’ = 3.9 (which means 
that the cell is young enough), and the concentration of the nitrogen 
compounds falls below the threshold value CT’ = 0.4. Production ps 
describes the subsequent growth of the heterocyst. 

Snapshots of the developmental sequence of Anabaena are given 
in Figure 1.35. The vegetative cells are shown as rectangles, colored 
according to the concentration of the nitrogen compounds (white means 
low concentration). The heterocysts are represented as red disks. The 
values of parameters C'H, CT and ST were selected to provide the 
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correct distribution of the heterocysts, and correspond closely to the 
values reported in [25]. The mathematical model made it possible to 
estimate these parameters, although they are not directly observable. 


1.10.3 Turtle interpretation of parametric words 


If one or more parameters are associated with a symbol interpreted by 
the turtle, the value of the first parameter controls the turtle’s state. 
If the symbol is not followed by any parameter, default values specified 
outside the L-system are used as in the non-parametric case. The basic 
set of symbols affected by the introduction of parameters is listed below. 


F(a) Move forward a step of length a > 0. The position of the 
turtle changes to (2’,y’, 2’), where x’ = 4+ aHz, y’ =y+aHy 
and 2’ = z+ad,. A line segment is drawn between points 
(x,y,z) and (2',y’, 2’). 


f(a) Move forward a step of length a without drawing a line. 


+(a) Rotate around U by an angle of a degrees. If a is positive, the 
turtle is turned to the left and if a is negative, the turn is to 
the right. 


& (a) Rotate around © by an angle of a degrees. If a is positive, 
the turtle is pitched down and if a is negative, the turtle is 
pitched up. 


/(a) Rotate around # by an angle of a degrees. If a is positive, the 
turtle is rolled to the right and if a is negative, it is rolled to 
the left. 


It should be noted that symbols +, &, A, and / are used both as 
letters of the alphabet V and as operators in logical and arithmetic 
expressions. Their meaning depends on the context. 

The following examples illustrate the operation of parametric L- 
systems with turtle interpretation. The first L-system is a coding of 
a Koch construction generating a variant of the snowflake curve (Fig- 
ure 1.1 on page 2). The initiator (production predecessor) is the hy- 
potenuse AB of a right triangle ABC (Figure 1.36). The first and the 
fourth edge of the generator subdivide AB into segments AD and DB, 
while the remaining two edges traverse the altitude CD in opposite di- 
rections. From elementary geometry it follows that the lengths of these 
segments satisfy the equations 


q=c-—p and h = ./pq. 


The edges of the generator can be associated with four triangles that 
are similar to ABC and tile it without gaps. According to the relation- 
ship between curve construction by edge rewriting and planar tilings 
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p D q 
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Figure 1.36: Construction of the generator for a “row of trees.” The edges 
are associated with triangles indicated by ticks. 


(Section 1.4.1), the generated curve will approximately fill the triangle 
ABC. The corresponding L-system is given below: 


#define c 1 

#define p 0.3 

#define q c—p 
#define h (pxq)A0.5 
ai (1) 


pi: F(a) F(a*p)+ F(w*h)——F(a*h) + F(w*q) 


The resulting curve is shown in Figure 1.37a. In order to better 
visualize its structure, the angle increment has been set to 86° instead 
of 90°. The curve fills different areas with unequal density. This results 
from the fact that all edges, whether long or short, are replaced by 
the generator in every derivation step. A modified curve that fills the 
underlying triangle in a more uniform way is shown in Figure 1.37b. It 
was obtained by delaying the rewriting of shorter segments with respect 
to the longer ones, as specified by the following L-system. 


w : F(1,0) 

pi: F(a,t):t=0 — F(x*p,2)+ F(xx*h,1)— 
—F(x*h,1)+ F(x *q,0) 

po: F(z,t):t>0 > F(z,t-1) 


The next example makes use of node rewriting (Section 1.4.2). The 
construction recursively subdivides a rectangular tile ABC'D into two 
tiles, AEF D and BCFE, similar to ABCD (Figure 1.38). The lengths 
of the edges form the proportion 
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Figure 1.37: Two curves suggesting a “row of trees.” Curve (b) is from [95, 


page 57]. 


which implies that b = a/\/2. Each tile is associated with a single-point 
frame lying in the tile center. The tiles are connected by a branching 
line specified by the following L-system: 


#define R 1.456 

w : A(1) (1.9) 

pi: A(s) > F(s)[+A(s/R)][-A(s/R)]| 
The ratio of branch sizes R slightly exceeds the theoretical value of 
/2. As a result, the branching structure shown in Figure 1.39 is self- 
avoiding. The angle increment was set arbitrarily to 6 = 85°. 

The L-system in equation (1.9) operates by appending segments 
of decreasing length to the structures obtained in previous derivation 
steps. Once a segment has been incorporated, its length does not 
change. A structure with identical proportions can be obtained by 
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Figure 1.38: Tiling associated with a space-filling branching pattern 





Figure 1.39: A branching pattern generated by the L-system specified in 
equation (1.9) 
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Figure 1.40: Initial sequences of figures generated by the L-systems specified 
in equations (1.9) and (1.10) 


appending segments of constant length and increasing the lengths of 
previously created segments by constant R in each derivation step. The 
corresponding L-system is given below. 


w: A 
pi: A — F(1)[+A][-A] (1.10) 
pz: F(s) — F(s*R) 


The initial sequence of structures obtained by both L-systems are 
compared in Figure 1.40. Sequence (a) emphasizes the fractal character 
of the resulting structure. Sequence (b) suggests the growth of a tree. 
The next two chapters show that this is not a mere coincidence, and 
the L-system specified in equation (1.10) is a simple, but in principle 
correct, developmental model of a sympodial branching pattern found 
in many herbaceous plants and trees. 


Chapter 2 


Modeling of trees 


Computer simulation of branching patterns has a relatively long history. 
The first model was proposed by Ulam [149], (see also [138, pages 127— 
131]), and employed the concept of cellular automata that had been 
developed earlier by von Neumann and Ulam [156]. The branching 
pattern is created iteratively, starting with a single colored cell in a 
triangular grid, then coloring cells that touch one and only one vertex 
of a cell colored in the previous iteration step. 

This basic idea gave rise to several extensions. Meinhardt [97, Chap- 
ter 15] substituted the triangular grid with a square one, and used the 
resulting cellular space to examine biological hypotheses related to the 
formation of net-like structures. In addition to pure branching patterns, 
his models capture the effect of branch reconnection or anastomosis that 
may take place between the veins of a leaf. Greene [54] extended cellu- 
lar automata to three dimensions, and applied the resulting voxel space 
automata to simulate growth processes that react to the environment. 
For instance, Figure 2.1 presents the growth of a vine over a house. 
Cohen [15] simulated the development of a branching pattern using ex- 
pansion rules that operate in a continuous “density field” rather than 
a discrete cellular or voxel space. 

The common feature of these approaches is the emphasis on inter- 
actions between various elements of a growing structure, as well as the 
structure and the environment. Although interactions clearly influence 
the development of real plants, they also add to the complexity of the 
models. This may explain why simpler models, ignoring even such fun- 
damental factors as collisions between branches, have been prevalent to 
date. The first model in that category was proposed by Honda [65] who 
studied the form of trees using the following assumptions (Figure 2.2). 


e Tree segments are straight and their girth is not considered. 


e A mother segment produces two daughter segments through one 
branching process. 
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Figure 2.1: Organic architecture by Greene [54] 


e The lengths of the two daughter segments are shortened by con- 
stant ratios, r; and ra, with respect to the mother segment. 


e The mother segment and its two daughter segments are contained 
in the same branch plane. The daughter segments form constant 
branching angles, a, and az, with respect to the mother branch. 


e The branch plane is fixed with respect to the direction of gravity 
so as to be closest to a horizontal plane.' An exception is made 
for branches attached to the main trunk. In this case, a constant 
divergence angle a between consecutively issued lateral segments 
is maintained (cf. Chapter 4). 


By changing numerical parameters, Honda obtained a wide vari- 
ety of tree-like shapes. With some improvements [38], his model was 
applied to investigate branching patterns of real trees [39, 66, 67, 68]. 
Subsequently, different rules for branching angles were proposed to cap- 
ture the structure of trees in which planes of successive bifurcations are 
perpendicular to each other [69]. The results of Honda served as a basis 
for the tree models proposed by Aono and Kunii [2]. They suggested 





‘More formally, the line perpendicular to the mother segment and lying in the 
branch plane is horizontal. 
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Figure 2.2: Specification of tree geometry according to Honda [65] 


several extensions to his model, the most important of which was the 
biasing of branch positions in a particular direction, applied to produce 
the effects of wind, phototropism and gravity. A similar concept was 
introduced previously by Cohen [15], while more accurate physically- 
based methods for branch bending were developed by de Reffye [28] 
and Armstrong [4]. 

The models of Honda and Aono and Kunii were rendered using 
straight lines of constant or varying width to represent “tree skele- 
tons.” A substantial improvement in the realism of synthetic images 
was achieved by Bloomenthal [11] and Oppenheimer [105], who intro- 
duced curved branches, carefully modeled surfaces around branching 
points, and textured bark and leaves (Figure 2.3). 

The approaches stemming from the work of Honda defined branching 
structures using deterministic algorithms. In contrast, stochastic mech- 
anisms are essential to the group of tree models proposed by Reeves and 
Blau [119], de Reffye et al. [30], and Remphrey, Neal and Steeves [120]. 
Although these models were described using different terminologies, 
they share the basic paradigm of specifying tree structures in terms of 
probabilities with which branches are formed. The work of Reeves and 
Blau aimed at producing tree-like shapes without delving into biological 
details of the modeled structures (Figure 2.4). In contrast, de Reffye 
et al. [29] used a stochastic approach to simulate the development of 
real plants by modeling the activity of buds at discrete time intervals. 
Given a clock signal, a bud can either: 


e do nothing, 
e become a flower, 


e become an internode terminated by a new straight apex and one 
or more lateral apices subtended by leaves, or 


e die and disappear. 
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Figure 2.4: A forest scene by Reeves [119] ©1984 Pixar 
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re: 


Figure 2.5: Oil palm tree canopy from CIRAD Modelisation Laboratory 


These events occur according to stochastic laws characteristic for each 
variety and each species. The geometric parameters, such as the length 
and diameter of an internode, as well as branching angles, are also 
calculated according to stochastic laws. 

The basic types of developmental rules incorporated into this method 
correspond to the 23 different types of tree architectures identified by 
Hallé, Oldeman and Tomlinson [58]. Detailed models of selected plant 
species were developed and are described in the literature [16, 20, 26, 
27, 76]. A sample tree model is shown in Figure 2.5. The approach of 
Remphrey [120, 121, 122] is similar to that of de Reffye, except that 
larger time steps are used (one year in the model of bearberry described 
in [120]). Consequently, the stochastic rules must describe the entire 
configuration of lateral shoots that can be formed over a one-year pe- 
riod. 

The application of L-systems to the generation of botanical trees was 
first considered by Aono and Kunii [2]. They referred to the original 
definition of L-systems [82] and found them unsuitable to model the 
complex branching patterns of higher plants. However, their arguments 
do not extend to parametric L-systems with turtle interpretation. For 
example, the L-system in Figure 2.6 implements those tree models of 
Honda [65] in which one of the branching angles is equal to 0, yielding 
a monopodial structure with clearly delineated main and lateral axes 
(see Section 3.2 for a formal characterization). 
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n = 10 

#define r; 0.9 /* contraction ratio for the trunk */ 
#define rp 0.6 /* contraction ratio for branches */ 
#define ap 45 /* branching angle from the trunk */ 
#define a) 45 /* branching angle for lateral axes */ 


#define d 137.5 /* divergence angle */ 
#define w, 0.707 /* width decrease rate */ 


w :A(1,10) 

pi: AC1l,w) :*— !(w)F(1) [& (ap) BC re, wew,) ]/(d) ACL ry , wew,) 
po: BCl,w) :*— !(w)F(1) [-(ag) $C (1*r2, wew,) ]C(1*ry, wwe) 
p3: C(l,w) :*— !(w)F(1) [+ (a2) $B(1*r2, wew,) ]B(1*ry , ww) 


Figure 2.6: Examples of the monopodial tree-like structures of Honda [65], 
generated using L-systems 
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According to production p;, the apex of the main axis A produces 
an internode F' and a lateral apex B in each derivation step. Con- 
stants r; and rg specify contraction ratios for the straight and lateral 
segments, ap and az are branching angles and d is the divergence angle. 
Module !(w) sets the line width to w, thus production p; decreases 
the width of the daughter segments with respect to the mother seg- 
ment by the factor w, = 0.707. This constant satisfies a postulate by 
Leonardo da Vinci [95, page 156], according to which “all the branches 
of a tree at every stage of its height when put together are equal in 
thickness to the trunk below them.” In the case where a mother branch 
of diameter w, gives rise to two daughter branches of equal diameter 
Wo, this postulate yields the equation w? = 2w3, which gives a value for 
w, equal to w2/w, = 1//2 ~ 0.707. A general discussion of the rela- 
tionships between the diameters of the mother and daughter branches 
is included in the book by Macdonald [94, pages 131-135]. 

Productions pz and p3 describe subsequent development of the lat- 
eral branches. In each derivation step, the straight apex (either B or C) 
issues a lateral apex of the next order at angle a2 or —dyg with respect 
to the mother axis. Two productions are employed to create lateral 
apices alternately to the left and right. The symbol $ rolls the turtle 
around its own axis so that vector © pointing to the left of the turtle 
(Section 1.5) is brought to a horizontal position. Consequently, the 
branch plane is “closest to a horizontal plane,” as required by Honda’s 
model. From a technical point of view, $ modifies the turtle orientation 
in space according to the formulae 


Vx 


Se and U=HxL, 
|V x A 


where vectors H, £ and U are the heading, left and up vectors associated 
with the turtle, V is the direction opposite to gravity, and |A| denotes 
the length of vector 4. The tree-like structures shown in Figure 2.6 
were generated using the values of constants listed in Table 2.1, and 
coincide with the structures presented by Honda. 








Figure | 7, | T2 | dao | ae 
a 0.9 | 0.6 | 45 | 45 
b 0.9|0.9 | 45 | 45 
Cc 0.9} 0.8 | 45 | 45 
d 0.9 | 0.7 | 30 | -30 























Table 2.1: Constants for the monopodial tree structures in Figure 2.6 
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A slightly different L-system, specified in Figure 2.7, captures sym- 
podial structures, where both daughter segments form non-zero angles 
with the mother segment. In this case the activity of the main apex is 
reduced to the formation of the trunk F' and a pair of lateral apices B 
(production p;). The subsequent branching pattern is captured by pro- 
duction pz. The sample structures in Figure 2.7 were obtained using the 
constants listed in Table 2.2, and correspond to the models presented 
by Aono and Kunii. 

The L-systems considered so far have been designed in a manner 
that emphasizes their relation to the models described in the litera- 
ture. Specifically, all segments are assigned their final length at the 
time of creation, and no further elongation occurs. As pointed out in 
Section 1.10.3, similar structures can be obtained by creating new seg- 
ments of constant length and increasing the lengths of previously cre- 
ated segments by a constant factor in each derivation step. A sample 
L-system constructed according to this paradigm is given in Figure 2.8. 

The overall structure of the tree is defined by production p,. In 
each derivation step, apex A produces three new branches terminated 
by their own apices. Parameter w and constant v, relate the width of 
the mother branch w, to that of the daughter branches w2. According 
to da Vinci’s postulate w? = 3w3, thus v, = w,/w2 = V3 & 1.732. 
Productions py and p3 capture the gradual elongation of branches and 
the increase in their diameter over time. 

The bending of branches is simulated by slightly rotating the tur- 
tle in the direction of a predefined tropism vector T after drawing each 
segment (Figure 2.9). The orientation adjustment a is calculated us- 
ing the formula a = e|H x T|, where e is a parameter capturing axis 
susceptibility to bending. This heuristic formula has a physical moti- 
vation; if T is interpreted as a force applied to the endpoint of vector #, 
and # can rotate around its starting point, the torque is equal to # x T. 
The parameters relevant to the generation of the tree-like structures in 
Figure 2.8 are listed in Table 2.3. A more realistic rendering of the tree 
in Figure 2.8d is presented in Figure 2.10. 








| Figure | ri] re | a | ae | 
a 0.9/0.7} 5 | 65 
b 0.9 | 0.7 | 10 | 60 
Cc 0.9 | 0.8 | 20 | 50 
d 0.9 | 0.8 | 35 | 35 

















Table 2.2: Constants for the sympodial tree structures in Figure 2.7 
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n = 10 


#define r, 0 /* contraction ratio 1 */ 
#define ro 0. /* contraction ratio 2 */ 
#define a, 10 /* branching angle 1 */ 
#define a 60 /* branching angle 2 */ 
#define w, 0.707 /* width decrease rate */ 


w : A(1,10) 

py : AC,w) : *—> !(w) FL) [&(a,) B(1 +r, , wtw,) J 
/ (180) [& (ap) B(1*r2, w*w,) J 

po : BCl,w) : *—!(w) FC) + (a,) $B(1*r,, wtw,)] 
[- (ap) $B(1*r2,wtw,) J 


Figure 2.7: Examples of the sympodial tree-like structures of Aono and 
Kunii [2], generated using L-systems 
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#define d, 94.74 /* divergence angle 1 */ 
#define d, 132.63 /* divergence angle 2 */ 
#define a 18.95 /* branching angle */ 
#define 1, 1.109 /* elongation rate */ 
#define v, 1.732 /* width increase rate */ 
w : !(1)F(200)/(45)A 


pi: A: * — !(vr)F (50) [&(a)F (50) Al / (dz) 
[& (a) F (50) A] / (do) [& (a) F (50) Al 

po : F(1) : *— F(1#1,) 

p3 : '!Cw) : *—> ! Cwrv,) 


Figure 2.8: Examples of tree-like structures with ternary branching 
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Figure 2.9: Correction a of segment orientation FH due to tropism T 





=> 


Figure dy dy a lL, T e€ 
94.74 | 132.63 | 18.95 | 1.109 | 0.00,-1.00,0.00 | 0.22 
137.50 | 137.50 | 18.95 | 1.109 | 0.00,-1.00,0.00 | 0.14 
112.50 | 157.50 | 22.50 | 1.790 | -0.02,-1.00,0.00 | 0.27 
180.00 | 252.00 | 36.00 | 1.070 | -0.61,0.77,-0.19 | 0.40 
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Table 2.3: Constants for the tree structures in Figure 2.8 





Figure 2.10: Medicine Lake by Musgrave et al. [101] 
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Figure 2.11: A surrealistic elevator 


Conclusions The examples given above demonstrate that the tree models of Honda, 
as well as their derivatives studied by Aono and Kunii, can be expressed 
using the formalism of L-systems. In a separate study, Shebell [130] also 
showed that L-systems can be applied to generate the architectural tree 
models of Hallé, Oldeman and Tomlinson [58]. These results indicate 
that L-systems may play an important role as a tool for biologically- 
correct simulation of tree development and synthesis of realistic tree 
images. However, the tree-like shapes created so far are rather generic 
(Figure 2.11), and models of particular tree species, directly based on 
biological data, are yet to be developed. L-systems have found more 
applications in the domain of realistic modeling of herbaceous plants, 
discussed in the next chapter. 


Chapter 3 


Developmental models of 
herbaceous plants 


The examples of trees presented in the previous chapter introduce L- 
systems as a plant modeling tool. They also illustrate one of the most 
striking features of the generative approach to modeling, called data 
base amplification [136]. This term refers to the generation of complex- 
looking objects from very concise descriptions — in our case, L-systems 
comprised of small numbers of productions. Yet in spite of the small 
size, the specification of L-systems is not a trivial task. 

In the case of highly self-similar structures, the synthesis methods 
based on edge rewriting and node rewriting are of assistance, as illus- 
trated by the examples considered in Section 1.10.3. However, a more 
general approach is needed to model the large variety of developmental 
patterns and structures found in nature. 

The methodology presented in this chapter is based on the simu- 
lation of the development of real plants. Thus, a particular form is 
modeled by capturing the essence of the developmental process that 
leads to this form. This approach has two distinctive features. 


e Emphasis on the space-time relation between plant parts. 
In many plants, organs in various stages of development can be 
observed at the same time. For example, some flowers may still 
be in the bud stage, others may be fully developed, and still 
others may have been transformed into fruits. If development 
is simulated down to the level of individual organs, such phase 
effects are reproduced in a natural way. 


e Inherent capability of growth simulation. The mathemati- 
cal model can be used to generate biologically correct images of 
plants at different ages and to create sequences of images illus- 
trating plant development in time. 
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The models are constructed under the assumption that organisms 
control the important aspects of their own development. According 
to Apter [3, page 44], this simplification must be accepted as a nec- 
essary evil, as long as the scope of the mathematical model is limited 
to an isolated plant. Consequently, this chapter focuses on the mod- 
eling and generation of growth sequences of herbaceous or non-woody 
plants, since internal control mechanisms play a predominant role in 
their development. In contrast, the form of woody plants is determined 
to a large extent by the environment, competition among branches and 
trees, and accidents [164]. 


3.1 Levels of model specification 


L-systems can be constructed with a variety of objectives in mind, rang- 
ing from a general classification of branching structures to detailed mod- 
els suitable for image synthesis purposes. Accordingly, the L-systems 
presented in this chapter are specified at three levels of detail. The most 
abstract level, called partial L-systems, employs the notation of nonde- 
terministic OL-systems to define the realm of possibilities within which 
structures of a given type may develop. Partial L-systems capture the 
main traits characterizing structural types, and provide a formal basis 
for their classification. Control mechanisms that resolve nondetermin- 
ism are introduced in the next level, termed L-system schemata.! The 
topology of individual plants and temporal aspects of their develop- 
ment are described at this level. Schemata are of particular interest 
from a biological point of view, as they provide an insight into the 
mechanisms that control plant development in nature. The geomet- 
ric aspects are added in complete L-systems that include information 
concerning growth rates of internodes, the values of branching angles, 
and the appearance of organs. The difference between all three levels is 
illustrated using models of a single-flower shoot as a running example. 


3.1.1 Partial L-systems 


Consider the development of a shoot which, after a period of vegetative 
growth, produces a single flower. The partial L-system is given below. 


oO. 2G 

pi: a> I1{Lla 
po: a > ILIA 
p3: AO K 


The lower-case symbol a represents the vegetative apex, while the upper- 
case A is the flowering apex, capable of forming reproductive organs. 


(3.1) 





‘Tn the literature, the term “scheme” is also used to denote the class of L-systems 
with the same alphabet and productions, but with different axioms [62, page 54]. 
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Figure 3.1: Single-flower shoot 


A derivation step corresponds to a plastochron, defined as the time 
interval between the production of successive internodes by the apex. 
At each step apex a has a choice of forming either leaf L, internode I 
and new apex a (production p,;), or forming the same structures and 
transforming itself into a flowering apex A (p2), which subsequently 
creates flower K (p3). Once this transformation or developmental switch 
has taken place it cannot be reversed, since there is no rule allowing 
the transformation of A to a. Examples of strings generated by the 
L-system specified in equation (3.1) are given below. 


I[L\A I[L]a I[Lla 
L\K  —soTLA[LJIA. CTL LJa 
L\K —oTLITILJK CTL TLIITLA 
TLL|K TULIP LI[LK 
T[L|I[LNI[LIK 


A diagrammatic representation of a single-flower inflorescence is shown 
in Figure 3.1. 


3.1.2 Control mechanisms in plants 


A partial L-system does not specify the moments in which develop- 
mental switches occur. The timing of these switches is specified at the 
level of L-system schemata, which incorporate mechanisms that con- 
trol plant development. In biology, these mechanisms are divided into 
two classes depending on the way information is transferred between 
modules. The term lineage (or cellular descent) refers to the transfer 
of information from an ancestor cell or module to its descendants. In 
contrast, interaction is the mechanism of information exchange between 
neighboring cells (for example, in the form of nutrients or hormones). 
Within the formalism of L-systems, lineage mechanisms are represented 
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by context-free productions found in OL-systems, while the simulation 
of interaction requires the use of context-sensitive 1L-systems and 2L- 
systems.” Several specific mechanisms are listed below. Although they 
are described from the modeling perspective, a relation to physiological 
processes observed in nature can often be found. 


Stochastic mechanism 


The simplest method for implementing a developmental switch is to 
use a stochastic L-system. In this case the vegetative apex a has a 
probability 7, of staying in the vegetative state, and 72 of transforming 
itself into a flowering apex A. 


WwW: a 

pi: aS I{Lla 
po: aS I[L|A 
p3 : WT 


The probability distribution (71, 72) is found experimentally, with 7+ 
Tg = 1. 


The effect of environment 


Many plants change from a vegetative to a flowering state in response 
to environmental factors such as temperature or the number of day- 
light hours. Such effects can be modeled using one set of productions 
(called a table) for some number of derivation steps, then replacing it 
by another set: 


Table 1 Table 2 
WwW: a pi: a >I [LA 
pi: a-I[Lja po: AvoK 


The concept of table L-systems (TOL-systems) was introduced and 
formalized by Rozenberg [62, 127]. Note that the use of tables provides 
only a partial solution to the problem of specifying the switching time, 
since a control mechanism external to the L-system is needed to select 
the appropriate table. 


Delay mechanism 


The delay mechanism operates under the assumption that the apex 
undergoes a series of state changes that postpone the switch until a 
particular state is reached. 





?The clarity of this dichotomy is somewhat obscured by parametric OL-systems, 
which can simulate the operation of context-sensitive L-systems using an infinite 
set of parameter values. 
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This is captured by the following L-system in the case of a single-flower 
shoot. 

Ww > ao 

Di : a; > I[Laiy1 O<i<n-1 

Pn : Qn 7 I[L|A 

Dnit: A wOK 


According to this model, the apex counts the leaves it produces. While 
it may seem strange that a plant would count, it is known that some 
plant species produce a fixed number of leaves before they start flow- 
ering. 


Accumulation of components 


A developmental mechanism based on the accumulation of components 
is similar to that of delay, but emphasizes the physiological nature of 
the counting process. According to this approach, counting is achieved 
by a monotonic increase or decrease in the concentration of certain cell 
components. This process can be captured by the following parametric 
L-system: 


w : a(0) 

pi: a(c) : c<C—-I[Lja(ct+ Ac) 
po: afc) : c2>C1[LlA 

p3: A ae ak —K 


(3.2) 


The parameter c indicates current concentration of the controlling com- 
ponents in the apex a. In each derivation step, this concentration is 
increased by a constant Ac. The developmental switch occurs when 
the concentration reaches the threshold value C’. 


Development controlled by a signal 


In many plants, the switch from a vegetative to a flowering state is 
caused by a flower-inducing signal transported from the basal leaves 
towards the apex. The time of signal initiation is determined using 
one of the previously described methods, for example by counting. A 
sample L-system is given below. 


wi: D(1)a(1) . . 

Pi: ali) : i<m—a(i+]) 
p2: a(i) t=m-— I[La(1) 
D3: Di): i<d — DUit+1) 
ve Di): t=d — S(1) 
Ds: S(t): i<u > S(t4+1) 
De: S(t): t=u re 

pr: S)<IF +: t=u > IS(1) 
pg: S(i)<a(j): * — I[L|A 
Dg: Ac* — K 


Crocus 
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The apex a produces internodes J and leaves L on the main axis (po). 
The time between the production of two consecutive segments, or the 
plastochron of the main axis, is equal to m derivation steps (p,). After 
a delay of d steps (p3), a signal S is sent from the plant base towards 
the apices (p4). This signal is transported along the main axis with 
a delay of uw steps per internode I (ps,p7). Production pg removes the 
signal from a node after it has been transported along the structure 
(€ stands for the empty string). When the signal reaches the apex, a is 
transformed into flowering state A (pg) which yields flower K’ (po). Note 
that the signal has to propagate faster than one node per plastochron 
(u < m), otherwise it would not be able to catch up with the apex. The 
above processes are illustrated by the following developmental sequence, 
ford = 4,m=2 andu=1. 


BOOS 


ae 
—~ 
ee 
—" 


Ba RS BS 


Although the above model may appear unnecessarily complicated, sig- 
nals are indispensable in the simulation of complex flowering sequences 
discussed later. 


3.1.3 Complete models 


The L-systems considered so far are not directly suitable for image syn- 
thesis purposes. To this end, they must be completed with geometric 
information. The relation between an L-system scheme and a corre- 
sponding complete L-system is discussed using the model of crocuses 
shown in Figure 3.2 as an example. 

The development is controlled using a delay expressed as an accu- 
mulation mechanism (equation (3.2)). In contrast to L-system schemes 
in which symbols represent module types, the L-system in Figure 3.2 
is specified in terms of turtle symbols. Production p,; describes the cre- 
ation of successive internodes F' and leaves L by the vegetative apex a. 
The leaves branch from the stem at an angle of 30° and spiral around 
the main axis with a divergence angle equal to 137.5° (see Chapter 4). 
Productions pz and p3 describe the developmental switch and the cre- 
ation of flower K taking place respectively in steps T,, and 7),,,. Pro- 
ductions p, and ps capture the development of leaves and flowers until 
they reach their final shapes 77, and 7, steps after creation. For each 
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#define T, 7 /* developmental switch time */ 
#define T, 9 /* leaf growth limit */ 
#define Tx 5 /* flower growth limit */ 
#include L(0),L(1),...,L(T,) /* leaf shapes */ 
#include K(0O),K(1),...,K(Tx) /* flower shapes */ 

a @ el) 

Pi : a(t) 2 t<T, — F(1) 18(30)<1L00)]/(137.5)a(t+1) 
po : a(t) : t=T, — F(20)A 

p3 : A : Ok — -«K(0) 

pa: L(t) : t<Tp — Ltt) 

ps : K(t) : t<Txy — K(t+1) 

pe : FC) : 1K2 - F(1+0.2) 


Figure 3.2: Crocuses 


Subapical 
growth 


Expression 
using 
L-systems 


70 Chapter 3. Developmental models 




















branch 
main apex main apex 
terminates continues 
all some all some 
lateral lateral lateral lateral 
apices apices apices apices 
terminate | continue terminate continue 
terminal | sympodial | monopodial | polypodial 














Table 3.1: Basic growth patterns of branching structures 


value of parameter t, the corresponding organ shapes L(t) and A(t) are 
modeled using bicubic patches incorporated into the plant structure as 
described in Section 5.1. Production pg specifies the gradual elongation 
of internodes. 


3.2. Branching patterns 


The partial L-system in equation (3.1) and the related schemata employ 
subapical growth mechanisms in which new branches are created exclu- 
sively by apices. All herbaceous plants develop this way. The archi- 
tecture of a branching structure is to a large extent determined by the 
relationships between terminal and continuing apices. While a contin- 
uing apex produces branches again and again, a terminal one either 
gives rise to an appendage such as a flower or dies. The possible com- 
binations are listed in Table 3.1. Of the four terms assigned to these 
possibilities, two are commonly used in biology, namely, sympodial and 
monopodial, while the other two terms are introduced here to denote 
the remaining cases, usually not characterized in the literature. 

The above branching patterns can be represented conveniently using 
partial bracketed L-systems. Let A,B,C denote continuing apices, X 
a terminal apex, and J an internode. The terminal and sympodial 
patterns are characterized by rules of the form 


Aa I[B]"[X]"X, 


with n = 0, m > 0 in the case of terminal patterns and n > 0, m > 0 
in the case of sympodial patterns. The important property is that the 
main apex terminates its development in all these cases. 

Monopodial and polypodial patterns have rules of the form 


A- I[B)"|X]"C, 
with n = 0, m > O in the case of monopodial patterns and n > 0, 


m > 0 in the case of polypodial ones. In these cases, it is important 
that the apex remains active. 
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The terms defined here apply to branching structures in general, 
whether they only produce vegetative organs (branches and leaves) or 
reproductive organs as well. In the latter case, the terminal organs 
develop from flower buds to flowers to fruits but do not give rise to 
vegetative structures (there are exceptions to this statement but they 
can be neglected here). 


3.3. Models of inflorescences 


The following discussion focuses on the modeling of compound flowering 
structures or inflorescences. In some cases an entire shoot system can 
be considered an inflorescence, in others only some of the branches 
bear flowers and are inflorescences. Inflorescence architecture is an 
elaboration of branching structures in general. 

In the domain of botanical applications of L-systems, the study 
of inflorescences has played a particularly visible role [44, 45, 46, 47, 
77, 86]. Unfortunately, the terms used for the various inflorescence 
types are not uniform in the literature. Besides a purely morphological 
terminology, attempts have been made to construct a “typological” 
terminology, expressing the “essential” features of flowering structures 
[144, 145, 157, 158]. However, these terms are not generally accepted 
[12]. A compromise has been proposed by D. and U. Miiller-Doblies 
[100], which serves as a basis for the classification that guides this 
presentation. 


3.3.1 Monopodial inflorescences 
Simple racemes (open) 


Racemes are characteristically monopodial inflorescences; a shoot has 
lateral apices with terminal structures and a main apex that continues 
to development. A raceme is open if the main apex does not form a 
flower. The partial L-system for this widely occurring type of inflores- 
cence is: 


case 
a eG TE (3.3) 


This system differs from that modeling a shoot with a single flower 
(equation (3.1)) only in production p3. Here it is designed to repeatedly 
produce lateral flowers (Figure 3.3), while in the previous system A 
produces a single flower. 


Classification 
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Figure 3.3: Open racemes: (a) elongated form, (b) planar form 





Figure 3.4: Lily-of-the-valley 
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The flowering sequence in open racemes is always acropetal (from 
base to top). This can be observed after substituting production p3 in 
the L-system in equation (3.3) with productions Ps and p4, which use 
indexed symbols K; to denote subsequent stages of flower development. 


pz: A > I[IKjJA 
pi hy => Tass 1>0 


The indexed notation K; — Kj, stands for a (potentially infinite) 
set of productions Ko — Ky, Ky — Ko, Ko — K3,.... The developmen- 
tal sequence begins as follows: 


Ty ILPKG]A 


At each developmental stage the inflorescence contains a sequence 
of flowers of different ages. The flowers newly created by the apex are 
delayed in their development with respect to the older ones situated at 
the stem base. Graphically, this effect is illustrated by the model of 
a lily-of-the-valley shown in Figure 3.4. The following quotation from 
d’Arcy Thompson [143] applies: 


A flowering spray of lily-of-the-valley exemplifies a growth- 
gradient, after a simple fashion of its own. Along the stalk 
the growth-rate falls away; the florets are of descending age, 
from flower to bud; their graded differences of age lead to 
an exquisite gradation of size and form; the time-interval 
between one and another, or the “space-time relation” be- 
tween them all, gives a peculiar quality — we may call it 
phase-beauty — to the whole. 


Another example of “phase beauty” can be seen in the shoot of shep- 
herd’s purse (Capsella bursa-pastoris) shown in Figure 3.5. Productions 
D1, Po and pz describe the activities of the apex in the vegetative and 
flowering states, in accordance with the L-system in equation (3.3). The 
developmental switch is implemented using a delay mechanism. Pro- 
ductions p, and p; capture the linear elongation of internodes in time, 
while pe and p7 describe the gradual increase of the angle at which 
the flower stalks branch from the main stem. Productions pg, po and 
pi specify the shapes of leaves L, flower petals K and fruits X using 
developmental surface models discussed in Section 5.2. Production pj 
controls the flowering time. Symbol % in the successor of production 
pi, simulates the fall of petals by cutting them off the structure at 
the time of fruit formation. The default value of the angle increment 
corresponding to the symbol + with no parameter is 18°. 


Acropetal 
sequence 


Lily-of-the- 
valley 


Capsella 
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Pi 
P2 
P3 


Pa 
Ps 
Pe 
P7 
Ps 


P9 


Pio + 
Pu: 


I(9)a(13) 

a(t) : t>0 — 
a(t) : t=0 — 
A * — 
I(t) : t>0 — 
I(t) : t=0 — 
u(t) : t>0 — 
u(t) : t=0 — 
L * — 
K * — 
X(t) : t>0 — 
X(t) : t=0 — 


[&(70)L] /(137.5)1(10)a(t-1) 
[&(70)L]/(137.5)1(10)A 
[&(18)u(4) FFI (10) 1(5)X (5) KKKK] 
/(137.5)1(8)A 

FI(t-1) 

F 

&(9)u(t-1) 

&(9) 

[{.-FI(7)+F1(7)+F1(7) }] 

[{ .+F1(7)-FI(7)-FI(7) }] 

[&{ .+F1I(2)--FI(2) }] 

[&{ .-FI(2)++F1(2) }J/(90) 
X(t-1) 

A (50) [[-GGGG++ [GGG[++G{.].].]. 
++GGGG.--GGG.--G.}1% 


Figure 3.5: Development of Capsella bursa-pastoris. Every fourth derivation 
step is shown. 
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Figure 3.6: Apple twig 


Simple raceme (closed) 


The inflorescence of an apple tree (Figure 3.6) provides an example of 
a closed raceme. In this case, the main apex eventually terminates its 
development and produces a terminal flower (Figure 3.7). The corre- 
sponding partial L-system is given below. 


Ww: 
pi: a > I[Lla 
po: a > ILIA 
p3: AvI|K]A 
pa: AWK 


Q 


Developmental switches are associated with two symbols, a and A. 
Thus, in order to obtain an L-system scheme it is necessary to specify 
how both of these switches will be controlled. 

The flowering sequence is usually acropetal but could also be basi- 
petal, z.e., progressing downward after the formation of the terminal 
flower on the main axis. In the latter case a basipetal signal, as dis- 
cussed in Section 1.8, can be applied to induce the transformation of 
dormant flower buds into flowers. 


Single-signal 
model 
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Figure 3.7: Closed racemes: (a) elongated form, (b) planar form 


Compound raceme (open dibotryoid) 


Racemes can also occur on complex branching structures. The sim- 
plest of these inflorescences is one with open racemes on the first order 
branches as well as on the main axis (Figure 3.8a). This two-level com- 
pound structure (thus dibotryoid) is described by the following partial 
L-system. 


Ww: 
Pi: 
P2: 
P3 : 
Pa : 
Ds : 
Deo : 
P7: |B 

Three developmental transformations are necessary: the first for the 
change from leaf to branch creation along the main axis (production po), 
the second for the switch from branching to lateral flower creation on 
the main axis (p4), and the third for the transition from leaf to lateral 
flower formation along the first-order branches (pg). Each branch is 
subtended by a leaf, which is why productions p3 and p, specify two 
appendages L and b. Branches with flowers K need not have subtending 
leaves, which is reflected in production p7. 

Within each component raceme, the flowering sequence is always 
acropetal, but the timing of switches has a crucial impact on the over- 
all flowering sequence and appearance of the plant. For example, let us 
assume that the switch from leaf to branch production is controlled by 
a delay, while the remaining two switches are caused by an acropetal 
flower-inducing signal (representing the hormone florigen). Such a de- 
velopment is captured by L-system 3.1 (see below). Initially, the veg- 
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a b 


Figure 3.8: Dibotryoids: (a) open, (b) closed 


#define d 13 /* delay for sending florigen */ 
#define e 3 /* delay for creating branches */ 
#define m 2 /* plastochron - main axis */ 
#define n 3 /* plastochron - lateral axis */ 
#define u 1 /* signal delay - main axis */ 
#define v 1 /* signal delay - lateral axis */ 
w :8(0)a(1,0) 

pi a(t,c) : t<m — a(tti,c) 

po: a(t,c) : (t=m)&(c<e) > 1(0,u) [L]a(1,c+1) 
Da a(t,c) : (t=m)&(c=e) > 1(0,u) [L] [b(1)]Ja(1,c) 
Haws BCG) ten — b(tt1) 

D5: b(t) :t=n — 1(0,v) [L]b(1) 

De: S(t) :* — S(tt+1) 

or? SG)-< TG, 4) stad — I11,j) 


pg: I(i,j) : (O<i)&(i<j) — IGi+1,j) 
po :1€i,j) < Ik,1) : G=j)&(k=0) — 14,1) 


Dio: I(i,j) < alk,1) : i>0 — I[L] [b(1)1B 
pu: I1(i,j) < b(K) :1>0 =. T [EIB 
Pi2: B 2k — I[K]B 


L-system 3.1: A model of dibotryoids 


Model analysis 
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etative apex a creates internodes J and leaves L with plastochron m 
(productions p; and pz). After the creation of e leaves a developmen- 
tal switch occurs, and apex a starts creating branches with the same 
plastochron (p3). The change of state is indicated by the value of the 
second parameter in the module a(t,c), which is now equal to e. The 
lateral apices b create internodes and leaves with plastochron n (p4 and 
ps). After a delay of d steps from the beginning of the simulation (pg), 
the flowering signal is introduced to the basal internode (p7), as indi- 
cated by a non-zero value of the first parameter in the module J (i, 7). 
The signal is passed along an axis at the rate of 7 steps per internode 
(pg and pg), where 7 = u for the main axis and j = v for the lateral 
axes. These rates are assigned to internodes by productions pz, p3 and 
ps. When the signal reaches an apex (either a or b), the apex is trans- 
formed into flowering state B (pio and pi;). From then on, new flowers 
K are produced in each derivation step (p12). 

In order to analyze the plant structure and flowering sequence re- 
sulting from the above development, let 7; denote the time at which 
apex b of the k-th lateral axis is transformed into the flowering state, 
and I; denote the length of this axis (expressed as the number of in- 
ternodes) at the transformation time. It is assumed here that the first e 
leaves count as lateral axes, thus k > e. Since it takes km time units to 
produce k internodes along the main axis and /,n time units to produce 
J, internodes on the lateral axis, we obtain: 


On the other hand, the transformation occurs when the signal reaches 
the apex. The signal is sent d time units after the development starts. 
It uses ku time units to travel through k zero-order internodes and /;,v 
time units to travel through /; first-order internodes: 


Ty =d+kut+tlkv 


Solving the above system of equations for J; and T), (and ignoring for 
simplicity some inaccuracy due to the fact that this system does not 
guarantee integer solutions), we obtain: 











Tei ey = 
n-v n—v 
pe, d 





n-v n-v 
In order to analyze the above solutions, let us first notice that the signal 
transportation delay v must be less than the plastochron of the lateral 
axes n, otherwise the signal would never reach the lateral apices. Under 
this assumption, the sign of the expression A = un—vm determines the 
overall flowering sequence, which is acropetal for A > 0 (Figure 3.9) 
and basipetal for A < 0 (Figure 3.10). If A = 0, all flowering switches 
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Figure 3.9: An acropetal flowering sequence in an open dibotryoid: m = 2, 
n=3,u=v=1, A=0.5; derivation lengths: 15—18—21—24—27—30-—33 


occur simultaneously. The sign of the expression m — wu determines 
whether the vegetative part of the shoot is more developed at the base 
(m —u <0) or near the top of the structure (m —u > 0). Figure 3.11 
shows a model of a member of the mint family that exhibits a basipetal 
flowering sequence. 


Compound racemes (closed dibotryoids) 


This inflorescence type differs from the previous one only in that each 
branch, including the main axis, bears a terminal flower (Figure 3.8b). 
A partial L-system can be obtained from that of equation (3.4) by 
adding one more production: 


p: BoK 
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Figure 3.10: A basipetal flowering sequence in an open dibotryoid: m = 2, 
n=5,u=1,v=3, A= —0.5; derivation lengths: 16 — 20 — 24 — 28 — 32 — 
36 — 40 
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Figure 3.11: A mint 


Compound raceme (closed tribotryoid) 


Racemic inflorescences can be compounded to a higher number of lev- 
els. The following is a partial L-system for a closed tribotryoid inflo- 
rescence, where closed racemes occur on second-order branches as well 
as on the terminal portions of first-order branches and of the main axis 
(Figure 3.12). The developmental process involves six developmental 


transformations. 


a ae pr: Bo IIM(4B 
po: a I|LiA ps : Bo I[L|[c|C 
psi Am I|L)[B)A a ae 
pa: A I|LI[O|B Oe Cea 
pe PUL eae 
pe: b> ILIB Digs Se 
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Figure 3.12: Closed tribotryoid 


3.3.2. Sympodial inflorescences 
Simple cymes (open) 


In racemes, the apex of the main axis produces lateral branches and 
continues to grow. In contrast, the apex of the main axis in cymes turns 
into a flower shortly after a few lateral branches have been initiated. 
Their apices turn into flowers as well, and second-order branches take 
over. In time, branches of higher and higher order are produced. Thus, 
the basic structure of a cymose inflorescence is captured by the partial 
L-system: 


1G 
pi: a-I[Lla 

po: av I{LiA (3.5) 
p3: Ao IlA|K 


As in the open raceme, there is a single symbol with alternative rules 
which specify that the vegetative apex a may change into a flower- 
producing apex A. Any one of the previously discussed mechanisms 
is available for timing this decision. Figure 3.13a shows an open cyme 
with branches curving in a spiral fashion, while Figure 3.13b shows one 
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b c 





Figure 3.13: Open cymes: (a) spiral form, (b) zig-zag form, (c) double 
with a zig-zag branching form. 


Double cymes (open) 


Frequently, not one but two lateral apices are produced under each 
terminal apex as in Figure 3.13c. In this case the partial L-system is: 


1 | ae I\Lja 
S a th (3.6) 
p3: Av I[AJ[Alk 


The two continuing lateral apices may develop at approximately 
equal rates (with the same plastochron) or with different rates, giving 
rise to asymmetric inflorescences. For example, the following L-system 
scheme describes the development of rose campion (Lychnis coronaria) 
as analyzed by Robinson [126]: 


WwW: Az 

Pr: Ay I[Ao] [Ay|l Ko 

p2: Aj; Ajat, 0<i<7 
Pat) KG Kia, 120 


Production p; shows that at their creation time, the lateral apices 
have different states Ag and Ay. Consequently, the first apex requires 
eight derivation steps to produce a flower and new branches, while the 
second requires only four steps. Each flower undergoes a sequence of 
changes, progressing from the bud stage to an open flower to a fruit. 
This developmental sequence is illustrated in Figure 3.14. According to 
production p,, the lateral apices branch at an angle of 45° and lie in a 
plane perpendicular to that defined by the mother axis and its sibling. 
Production p3 describes the linear elongation of internodes, while p4 


Lychnis 
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#include L(0),L(1),... /* leaf shapes */ 

#include K(0),K(1),... /* flower shapes */ 

w : A(T) 

pi: A(t) : t=7 — FI(20) [& (60) ~L(0) J / (90) [& (45) ACO) J / (90) 
[& (60) ~L(0)]/ (90) [&(45) A(4) ] FI (10) ~K (0) 

po: A(t) : t<7 > A(t+1) 

p3: I(t) : t>0 > FFI(t-1) 

pa: L(t): * Ltt) 

ps: K(t):* —K(tt1) 


Figure 3.14: Development of Lychnis coronaria 
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a b Cc 


Figure 3.15: Thyrsus: (a) spiral form, (b) zig-zag form, (c) double 


and ps capture the development of leaves and flowers over time. It is 
interesting to note that at different developmental stages there are some 
open flowers that have a relatively uniform distribution over the entire 
plant structure. This is advantageous to the plant since it increases the 
time span over which seeds will be produced. 


Cymes (closed) 


Sympodial inflorescences that produce a terminal flower at some point 
during their development are called closed cymes. They result from the 
addition of production 

pa: Aw Kk 


to the L-systems specified in (3.5) and (3.6), which define open single 
and open double cymes. 


Thyrsus (closed) 


A thyrsus is an inflorescence with branches of cymes borne on a monopo- 
dially branching axis. Thus, it represents a mixed sympodial and 
monopodial organization. Depending on the orientation of the flow- 
ers, a distinction between a thyrsus with cymes in a spiral form and in 
a zig-zag form can be made (Figure 3.15, a and b). Both of these types 
are described by the following partial L-system: 
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BSG 
pi: aol[Lia 
po: al |LiA 


p3: A= IIL|[BJA 


pa: Ao Kk 
ps: BoIlBlK 
pep: Bo K 


In addition, a thyrsus may have double cymes (Figure 3.15c). In the 
closed structure there are three developmental transformations. The 
first represents the change from vegetative to flowering development on 
the main axis (production pz). The second is necessary for the closure 
of the main axis with a terminal flower (p,). Both switches are related 
to the monopodial development of the main axis. The third transfor- 
mation is responsible for the formation of the flowers that terminate 
the development of the sympodial structures (pg). 


3.3.3 Polypodial inflorescences 
Panicle 


The term polypodial is not used in the botanical literature but is coined 
here to draw attention to the type of branching that represents con- 
tinuing development of the main axis as well as of the lateral apices 
of a branch. The corresponding inflorescence type is usually called a 
panicle. The presence of two continuing apices at each new node is 
expressed by the following production: 


A— I[L][AJA 


Since there can be nodes near the base of the plant that do not bear 
branches, the usual initial rules are included to model the transition 
from a purely vegetative to a flowering state. The resulting partial 
L-system is: 


OG 
pi: a—I[Lla 
po: a I[L\A 
p3: A I|L|[A]A 
pa: AO K 


An example of a paniculate structure is shown in Figure 3.16. Note 
the presence of higher order branching and the lack of terminal racemes. 
Due to the repetitive application of production p3 at various levels of 
branching, the resulting structure is highly self-similar. The model in- 
cludes only two types of developmental transformations: the switch 
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Figure 3.16: Panicle (elongated form) 


from purely vegetative growth to the formation of the branching struc- 
ture (production pz), and the creation of terminal flowers (p4). The 
timing of the last production determines the flowering sequence of the 
plant. Two possible control mechanisms will be examined in detail, us- 
ing developmental models of the branching part of wall lettuce (Mycelis 
muralis) as examples. 

The development of Mycelis is difficult to model for two reasons. 
First, the plant exhibits a basipetal flowering sequence, which means 
that flowering starts at the top of the plant and proceeds downwards. 
Secondly, at some developmental stages the plant has an acrotonic 
structure, where the upper branches are more developed than the lower 
ones. Both phenomena are in a sense counter-intuitive, since it would 
seem that the older branches situated near the plant base should start 
growing and producing flowers before the younger ones at the plant 
top. To explain these effects, several models were proposed and for- 
mally analyzed by Janssen and Lindenmayer [77]. Their model IJ is 
restated here as parametric L-system 3.2. 

The axiom consists of three components. Modules F’ and A(0) rep- 
resent the initial segment and the apex of the main axis. Module /(20) 
is the source of a signal representing florigen. In nature, florigen is 
sent towards the apex by leaves located at the plant base, which is not 
included in this model. 

The developmental process consists of two phases that take place 
along the main axis and are repeated recursively in branches of higher 
orders. First, the main axis is formed in a process of subapical growth 
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#include 0 /* flower shape specification */ 
#ignore / + ~ 0 


w : I1(20)FA(0) 


Pat S < A(t) >* — T(O)<0 

po: A(t) :t>0— A(t-1) 

D3: A(t) :t=0— [+(30)G] F/(180) A(2) 
pa: S <F :* > FS 

D5: F>T(c) :* —- T(ct1)FU(c-1) 

pe : Ulc) <G :* —> I(c)FA(2) 

Dr: TCG) :c>0— I(c-1) 

Dg : I(c) :c=0— § 

po : S DK > € 

P10: Tee) :*# > € 


L-system 3.2: Mycelis muralis — Model II 


specified by production p3. The apex produces consecutive segments 
F at the rate of one segment every three derivation steps (the delay is 
controlled by production pz), and initiates branches G positioned at an 
angle of 30° with respect to the main axis. The symbol G is interpreted 
here in the same way as F’. At this stage, the branches do not develop 
further, which simulates the effect of apical dominance or the inhibition 
of branch development during the active production of new branches 
by the apex. 

After a delay of 20 derivation steps, counted using production pz, 
an acropetal flower-inducing signal S is sent by production pg. Produc- 
tion p4 transports S across the segments at the rate of one internode 
per step. Since new internodes are produced by the apex at a three 
times slower rate, the signal eventually reaches the apex. At this point, 
the second developmental phase begins. Production p,; transforms apex 
A(t) into a bud O. Further branch production is stopped and a signal 
T(c) is sent towards the base in order to enable the development of 
lateral branches. Parameter c is incremented by production ps each 
time signal T(c) traverses an internode. Subsequently, production pg 
introduces the value of parameter c into the corresponding branches, 
using module U(c) as a carrier. The successor of production pg has 
the same format as the axiom, thus module J(c) determines the delay 
between the initiation of branch development and time signal 5S, sent 
to terminate further internode creation. This delay c is smallest for 
the top branches and increases towards the plant base. Consequently, 
parameter c can be interpreted as the growth potential of the branches, 
allowing lower branches to grow longer than the higher ones. On the 
other hand, the development of the upper branches starts sooner, thus 
in some stages they will be more developed than the lower ones, and 
the flowering sequence will progress downwards, corresponding to ob- 
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#include K 
#consider MS T V 


™~ 


* flower shape specification */ 


w : 1(20)FA(0) 


pi: S < A(t) :* — TV K 
po : V< A(t) :* — TV K 
D3: A(t) :t>0 — A(t-1) 
Da: A(t) :t=0 — M[+(30)G]F/(180)A(2) 
pp: S <M 1k — 8 

pe : S > T:* — T 

pr: T<G 1 — FA(2) 
pe: V<M 1k — 8 

pg : T>V:* — W 

Pio: W ot =r ¥ 

Pi: I(t) t>0 — I(t-1) 
P12: I(t) :t=-0 —S§ 


L-system 3.3: Mycelis muralis — Model III 


servations of the real plant [77]. 

A diagrammatic developmental sequence of Mycelis muralis simu- 
lated using L-system 3.2 is shown in Figure 3.17. Initially, the segments 
are shown as bright green. The passage of florigen S turns them purple, 
and the lifting of apical dominance changes their color to dark green. 
Figure 3.18 represents a three-dimensional rendering of the same model. 
The three-dimensional structure differs from the two-dimensional di- 
agram only in details. The angle value associated with the module 
“/” in production p3 has been changed to 137.5°, resulting in a spiral 
arrangement of lateral branches around the mother axis. The leaves 
subtending branches have been included in the model, and flowers have 
been assumed to undergo a series of changes from bud to open flower 
to fruit. 

Another developmental model of Mycelis, referred to here as model 
III, is given by L-system 3.3. The initial phases of development are the 
same as in model II. First, apex A creates the main axis and initiates 
lateral branches (productions p3 and py). Symbol M in the successor of 
production p, marks consecutive branching points. After a delay of 20 
steps (w) counted by production p,1, flowering signal S is generated at 
the inflorescence base (pi2) and sent up along the main axis (ps). Upon 
reaching the apex, S induces its transformation into a terminal flower 
K, and initiates two basipetal signals T and V (p,). The basipetal 
signals also can be initiated by production py, which is needed for the 
proper timing of signals in the topmost lateral branch. Signal T’ prop- 
agates basipetally at the rate of one internode per derivation step (p¢) 
and lifts apical dominance, thus allowing the lateral branches to grow 
(p7). The presence of the second basipetal signal V is the distinctive 


Model III 
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Figure 3.17: Development of Mycelis muralis 
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Figure 3.18: A three-dimensional rendering of the Mycelis model 


feature of model III. Its role is to enable the formation of flowers on 
the lateral branches by generating the flowering signal S at their bases 
(pg). Since signal V propagates down the main axis at the rate of one 
internode per two derivation steps (p9, Pio), the interval between the 
lifting of apical dominance by signal 7 and induction of flowering signal 
S by signal V increases linearly towards the inflorescence base. This 
allows the lower branches to grow longer than the upper ones, resulting 
in a structure that is more developed near the base than near the apex 
in later developmental stages. 

This entire control process repeats recursively for each axis: its apex 
is transformed into a flower by signal S, the growth of lateral axes is 
successively enabled by signal 7’, and the second basipetal signal V is 
sent to induce the flowering signal S in the next-order axes. Conse- 
quently, a basipetal flowering sequence is observed along all axes of the 
panicle. 

Model II controls the flowering on lateral branches using growth 
potential c accumulated by signal 7’ on its way down, while model 
III employs the time interval between signals T’ and V for the same 
purpose. Since both models produce identical developmental sequences, 
it is not possible to decide which one is more faithful to nature without 
gathering additional data related to plant physiology. Nevertheless, the 
models clearly indicate that the flowering sequence of Mycelis cannot 
be explained simply in terms of two commonly recognized mechanisms, 
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n=10, 6 = 60° 


#include K /* flower shape specification */ 


w : AnvxK 


po: (OA : Ox — [-/~K] [+/~K]1(0)/(90)A 
po: Itt) : !(t=2) — FI(t+1) 
pst: Te) 2 G=2 — I(t+1) [-FFA] [+FFA] 


Figure 3.19: Lilac inflorescences 
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Figure 3.20: Geometry of a lilac inflorescence: (a) the decussate branching 
pattern, (b) the inflorescence skeleton without flowers 


the flowering signal and the lifting of apical dominance. Another factor, 
whether it is an accumulated delay or a third signal, is needed. The 
mathematical models bring forward evidence and assist in formulating 
plausible hypotheses related to the control mechanisms that may be 
employed by nature. The final answer will require further study of the 
real plant. 

The models of Mycelis employ relatively complicated control pro- 
cesses to explain the developmental sequence of a plant. On the other 
hand, if only a static image of a panicle in a particular developmental 
stage is needed, much simpler L-systems can be used. The L-system 
that generates the lilac inflorescences shown in Figure 3.19 is an ex- 
ample. Production p; describes the subapical development of an axis. 
Production pz models linear elongation of internodes in time and in- 
troduces a delay before p3 creates the lateral axes. The rotation of the 
apex by 90° (p,) results in a decussate branching pattern with consec- 
utive pairs of (n + 1)-order axes lying in the planes that pass through 
the n-order axis and are perpendicular to each other (Figure 3.20). A 
scene including lilac inflorescences is shown in Figure 3.21. 


3.3.4 Modified racemes 


There are four frequently encountered types of inflorescences that are 
morphological modifications of racemes. Their mature forms are of a 
special kind and need to be specified separately. 


Lilac 
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Figure 3.21: The garden of L 
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Figure 3.22: Umbels: (a) simple, (b) compound 


Umbel 


An umbel is characterized by more than two internodes attached to a 
single node, resulting in a typical umbrella-like shape. In a simple um- 
bel there are flowers at the ends of the lateral internodes (Figure 3.22a), 
while in compound umbels the branching pattern is repeated recursively 
a certain number of times (Figure 3.22b). The partial L-system for a 
simple umbel is 

w: A 

pp: Aw ILK]" 


and for a compound umbel of recursion depth two is 


w: A 

pi: A> I[IB|FB 
po: B= I{IC}!C 
p3: C—I[lk]™ 


This type of inflorescence is commonly found in the family Umbelliferae. 
For example, Figure 3.23 presents a model of a wild carrot. Note that 
the size of leaves decreases towards the top of the plant, producing a 
phase effect similar to that observed in simple racemes. In contrast, the 
most developed inflorescence is placed at the top of the plant, indicating 
developmental control by a hormone similar to that observed in mints 
(Figure 3.11). 


Spike 


An elongated raceme with closely packed flowers is called a spike. Many 
grasses and sedges have this kind of inflorescence (Figure 3.24a). See 
Figure 4.17 (page 117) for a realistic model. 


Wild carrot 


96 Chapter 3. Developmental models 





Figure 3.23: Wild carrot 


Spadix 


A fleshy elongated raceme is called a spadiz, and is frequently found in 
the family Araceae (Figure 3.24b). 


Capitulum 


A fleshy spherical or disk-shaped raceme is called a capitulum or “head.” 
The sunflower head is an inflorescence of this kind, the oldest flowers be- 
ing at the margin and the youngest at the center (Figure 3.24c). Mem- 
bers of the family Compositae commonly have this type of structure. 
One characteristic feature is the spatial arrangement of components, 
such as flowers or seeds, which form early discernible spiral patterns. A 
detailed description of these patterns is presented in the next chapter. 
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Figure 3.24: Modified racemes: (a) spike, (b) spadix, (c) capitulum 


Chapter 4 


Phyllotaxis 


The regular arrangement of lateral organs (leaves on a stem, scales on 
a cone axis, florets in a composite flower head) is an important aspect 
of plant form, known as phyllotaxis. The extensive literature generated 
by biologists’ and mathematicians’ interest in phyllotaxis is reviewed 
by Erickson [36] and Jean [78]. The proposed models range widely from 
purely geometric descriptions (for example, Coxeter [17]) to complex 
physiological hypotheses tested by computer simulations (Hellendoorn 
and Lindenmayer [59], Veen and Lindenmayer [151], Young [163]). This 
chapter presents two models suitable for the synthesis of realistic images 
of flowers and fruits that exhibit spiral phyllotactic patterns. 

Both models relate phyllotaxis to packing problems. The first oper- 
ates in a plane and was originally proposed by Vogel [154] to describe 
the structure of a sunflower head. A further detailed analysis was given 
by Ridley [124, 125]. The second model reduces phyllotaxis to the 
problem of packing circles on the surface of a cylinder. Its analysis was 
presented by van Iterson [75] and reviewed extensively by Erickson [36]. 

The area of phyllotaxis is dominated by intriguing mathematical 
relationships. One of them is the “remarkable fact that the numbers of 
spirals which can be traced through a phyllotactic pattern are predom- 
inantly integers of the Fibonacci sequence” [36, p. 54]. For example, 
Coxeter [17] notes that the pineapple displays eight rows of scales slop- 
ing to the left and thirteen rows sloping to the right. Furthermore, it is 
known that the ratios of consecutive Fibonacci numbers F;,.41/F, con- 
verge towards the golden mean Tt = (5+ 1)/2. The Fibonacci angle 
360°7~?, approximately equal to 137.5°, is the key to the first model 
discussed in this chapter. 


Vogel’s formula 
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#define a 137.5 /* divergence angle */ 
#include D /* disk shape specification */ 


w : ACO) 
py : Atm) : * > +(a) [£(mA0.5)~D] ACm+1) 


Figure 4.1: Pattern of florets in a sunflower head, according to Vogel’s for- 
mula 


4.1 The planar model 


In order to describe the pattern of florets (or seeds) in a sunflower head, 
Vogel [154] proposed the formula 


@=n*1387.5°, r=c/n, (4.1) 


where: 


e n is the ordering number of a floret, counting outward from the 
center. This is the reverse of floret age in a real plant. 


e @ is the angle between a reference direction and the position vec- 
tor of the n‘” floret in a polar coordinate system originating at 
the center of the capitulum. It follows that the divergence an- 
gle between the position vectors of any two successive florets is 
constant, a = 137.5°. 


e r is the distance between the center of the capitulum and the 
center of the n“” floret, given a constant scaling parameter c. 
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Figure 4.2: Generating phyllotactic patterns on a disk. These three patterns 
differ only by the value of the divergence angle a, equal to (a) 137.3°, (b) 
137.5° (the correct value), and (c) 137.6°. 


The distribution of florets described by formula (4.1) is shown in Fig- 
ure 4.1. The square-root relationship between the distance r and the 
floret ordering number n has a simple geometric explanation. Assum- 
ing that all florets have the same size and are densely packed, the total 
number of florets that fit inside a disc of radius r is proportional to 
the disk area. Thus, the ordering number n of the most outwardly 
positioned floret in the capitulum is proportional to r?, or r ~ \/n. 

The divergence angle of 137.5° is much more difficult to explain. 
Vogel [154] derives it using two assumptions. 


e Each new floret is issued at a fixed angle a with respect to the 
preceding floret. 


e The position vector of each new floret fits into the largest existing 

gap between the position vectors of the older florets. 

Ridley [125] does not object to these basic assumptions, but indi- 
cates that they are insufficient to explain the origin of the Fibonacci 
angle, and points to several arbitrary steps present in Vogel’s deriva- 
tion. He describes the main weakness as follows: 


While it is reasonable to assume that the plant could con- 
tain genetic information determining the divergence angle 
to some extent, it is completely impossible for this alone to 
fix the divergence angle to the incredible accuracy occurring 
in nature, since natural variation in biological phenomena 
is normally rather wide. For example, for the 55- and 89- 
parastichies to be conspicuous, as occurs in most sunflower 
heads, d must lie between = and x a relative accuracy of 
one part in 1869. 


The critical role of the divergence angle a is illustrated in Figure 4.2. 


Model 
justification 


Parastichies 
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Figure 4.3: Close-up of a daisy capitulum 


Although a comprehensive justification of Vogel’s formula may require 
further research, the model correctly describes the arrangement of flo- 
rets visible in actual capitula. The most prominent feature is two sets 
of spirals or parastichies, one turning clockwise, the other counterclock- 
wise, which are composed of nearest neighboring florets. The number 
of spirals in each set is always a member of the Fibonacci sequence; 21 
and 34 for a small capitulum, up to 89 and 144 or even 144 and 233 for 
large ones. For example, the capitulum of a daisy (Figure 4.3) exhibits 
34 clockwise spirals and 21 counterclockwise spirals, with directions 
determined by following the spirals outwards from the capitulum cen- 
ter. In the image of a domestic sunflower capitulum (Figure 4.4), one 
can discern 34 spirals running clockwise and 55 spirals running counter- 
clockwise. The number of perceived spirals depends on the capitulum 
size expressed in terms of the number of component florets. If the field 
of attention is limited to a circle approximately 2/3 the size of the en- 
tire sunflower capitulum in Figure 4.4, the number of discernible spirals 
becomes 34 and 21. 


Figure 4.4: Domestic sunflower head — 


Al. 


The planar model 


#define S 
#define R 
Ww A(O) 
py : AC) 
poy : Cn) 
p3 : C(n) 
pe ¢ CO) 
ps : C(n) 
pe : Cn) 
pz : C(n) 


/* seed shape */ 
/* ray floret shape */ 
#include MN OP /x petal shapes */ 


* — +(137.5) [£ (nA0.5)C(n)J]A(n+1) 
n <= 440 
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Sunflower head 
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The dependence of the number of parastichies on the size of the 
field of attention is yet another intriguing aspect of spiral phyllotaxis, 
as pointed out in the following excerpt from a letter by Alan Turing ! 
quoted in [18]: 


According to the theory I am working on now there is a 
continuous advance from one pair of parastichy numbers to 
another, during the growth of a single plant. ... You will 
be inclined to ask how one can move continuously from one 
integer to another. The reason is this — on any specimen 
there are different ways in which the parastichy numbers can 
be reckoned; some are more natural than others. During the 
growth of a plant the various parastichy numbers come into 
prominence at different stages. One can also observe the 
phenomenon in space (instead of in time) on a sunflower. 
It is natural to count the outermost florets as say 21 + 34, 
but the inner ones might be counted as 8+ 13. ... I don’t 
know any really satisfactory account, though I hope to get 
one myself in about a year’s time. 


A complete model of a flower head, suitable for realistic image syn- 
thesis, should contain several organs of various shapes. This is easily 
achieved by associating different surfaces with specific ranges of the in- 
dex n. For example, consider the L-system that generates the sunflower 
head (Figure 4.4). The layout of components is specified by produc- 
tion p;, similar to that of the L-system in Figure 4.1. Productions po 
to pz determine colors and shapes of components as a function of the 
derivation step number. The entire structure shown in Figure 4.4 was 
generated in 630 steps. Alternatively, random selection of similar sur- 
faces could have been employed to prevent the excessive regularity of 
the resulting image. 

Other extensions to the basic model consist of varying organ orien- 
tation in space and changing their altitude from the plane of the head 
as a function of n. For example, the sunflower plants included in Fig- 
ure 4.5 have flowers in four developmental stages: buds, young flowers 
starting to open, open flowers and older flowers where the petals begin 
to droop. All flowers are generated using approximately the same num- 
ber of florets. The central florets are represented by the same surface at 
each stage. The shape and orientation of surfaces representing petals 
vary from one stage to another. The plants have been modeled as di- 
botryoids, with a single signal inducing a basipetal flowering sequence, 
as described in the previous chapter. 





'To computer scientists, Alan Turing is best known as the inventor of the Turing 
machine [146], which plays an essential role in defining the notion of an algorithm. 
However, biologists associate Turing’s name primarily with his 1952 paper, “The 
chemical basis of morphogenesis” [147], which pioneered the use of mathematical 
models in the study of pattern formation and advocated the application of comput- 
ers to simulate biological phenomena. 
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Figure 4.5: Sunflower field 


Other examples 
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Figure 4.6: Zinnias 


The zinnias (Figures 4.6 and 4.7) illustrate the effect of changing a 
petal’s altitude, size and orientation as a function of n. The height at 
which a petal is placed decreases by a small amount as n increases. The 
size of each successive petal is incremented linearly. The orientation is 
also adjusted linearly by a small angle increment. Thus, petals with 
small values of index n are placed more vertically, while petals with 
larger indices n are more horizontal. Although the family Compositae 
offers the most examples of phyllotactic patterns, the same model can 
be applied to synthesize images of other flowers, such as water-lilies 
(Figures 4.8 and 4.9) and roses (Figure 4.10). 


Figure 4.7: Close-up of zinnias —> 
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Figure 4.8: Water-lily 








Figure 4.9: Lily pond 
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Figure 4.10: Roses 


4.2 The cylindrical model 


The spiral patterns evident in elongated organs such as pine cones, 
fir cones and pineapples, can be described by models that position 
components, in this case scales, on the surface of a cylinder. Van Iterson 
[75] divides phyllotactic patterns on cylinders into simple and conjugate 
ones. In the case of a simple arrangement, all components lie on a single 
generative helix. In contrast, conjugate patterns consist of two or more 
interleaved helices. This paper discusses simple phyllotactic patterns 
only. They are generally characterized by the formula 


C=2#FO, r= const, H=h*¥n, (4.2) 
where: 


e nis the ordering number of a scale, counting from the bottom of 
the cylinder; 


e ¢, r and H are the cylindrical coordinates of the n‘” scale; 


e a is the divergence angle between two consecutive scales (as in 
the planar case, it is assumed to be constant); and 


e h is the vertical distance between two consecutive scales (mea- 
sured along the main axis of the cylinder). 


Basic model 


Implementation 
using 
L-systems 


Analysis of 
model geometry 
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A parametric L-system that generates the pattern described by for- 
mula (4.2) is given in Figure 4.11. The operation of this L-system 
simulates the natural process of subapical development characterized 
by sequential production of consecutive modules by the top part of the 
growing plant or organ. The apex A produces internodes f(h) along 
the main axis of the modeled structure. Associated with each intern- 
ode is a disk ~D placed at a distance r from the axis. This offset is 
achieved by moving the disk away from the axis using the module f(r), 
positioned at a right angle with respect to the axis by &(90). The spi- 
ral disk distribution is due to the module /(a), which rotates the apex 
around its own axis by the divergence angle in each derivation step. 

In the planar model, the constant divergence angle a = 137.5° is 
found across a large variety of flower heads. The number of perceived 
parastichies is determined by the capitulum size, and it changes as the 
distance from the capitulum center increases. In contrast, a phyllotactic 
pattern on the surface of a cylinder is uniform along the entire cylinder 
length. The number of evident parastichies depends on the values of 
parameters a and h. The key problem, both from the viewpoint of 
understanding the geometry of the pattern and applying it to generate 
synthetic images, is to express the divergence angle a and the vertical 
displacement h as a function of the numbers of evident parastichies en- 
circling the cylinder in the clockwise and counterclockwise directions. A 
solution to this problem was proposed by van Iterson [75] and reviewed 
by Erickson [36]. Our presentation closely follows that of Erickson. 

The phyllotactic pattern can be explained in terms of circles packed 
on the surface of the cylinder. An evident parastichy consists of a 
sequence of tangent circles, the ordering numbers of which form an 
arithmetic sequence with difference m. The number m is referred to as 
the parastichy order. Thus, the circles on the cylinder surface may be 
arranged in two congruent 2-parastichies, five congruent 5-parastichies, 
and so on. The angular displacement between two consecutive circles 
in an m-parastichy is denoted by 6,,. By definition, 6,, belongs to the 
range (—7,7] radians. The relation between the angular displacement 
dm and the divergence angle a is expressed by the equation 


dm = ma — A,p27, (4.3) 


where A,, is an integer called the encyclic number. It is the number 
of turns around the cylinder, rounded to the nearest integer, which 
the generative helix describes between two consecutive points of the 
m-parastichy. 

Usually, one can perceive two series of parastichies running in oppo- 
site directions (Figure 4.11). The second parastichy satisfies an equa- 
tion analogous to (4.3): 


On = na — A,2r (4.4) 


Consider the m- and n-parastichies starting at circle 0. In their paths 
across the cylinder, they will intersect again at circle mn. Assume 
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#define a 137.5281 /* divergence angle */ 


#define h 35.3 /* vertical displacement */ 

#define r 500 /* component offset from main axis */ 
#include D /* disk shape specification */ 

wi: A 


py: A: * > [&(90)f(r)~D]f(h)/(a)A 


Figure 4.11: Parastichies on the surface of a cylinder and on the unrolled 
cylinder. The L-system generates the cylindrical pattern. 


that m and n are relatively prime; otherwise the phyllotactic pattern 
would have to contain several circles lying at the same height H and, 
contrary to the initial assumption, would not be simple. The circle 
myn is the first point of intersection between the m-parastichy and the 
n-parastichy above circle 0. Consequently, the path from circle 0 to mn 
along the m-parastichy, and back to 0 along the n-parastichy, encircles 
the cylinder exactly once. The section of m-parastichy between circles 0 
and mn consists of n+1 circles (including the endpoints), so the angular 
distance between the circles 0 and mn is equal to nd,,. Similarly, the 
distance between circles 0 and mn, measured along the n-parastichy, 
can be expressed as md,. As a result, 


Nbm — Moy = £2. (4.5) 
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Figure 4.12: An opposite parastichy triangle (as in Erickson [36, Fig. 3.8]). 
The base is formed by the circumference of the cylinder. The sides are 
formed by the parastichies. 


The signs in equation (4.5) correspond to the assumption that the spi- 
rals encircle the cylinder in opposite directions; thus one of the values 
6 is positive and the other one is negative. Substituting the right sides 
of equations (4.3) and (4.4) for 6,, and 6, yields 


nA, — mA, = +1. (4.6) 


To further analyze the pertinent geometric relationships, the cylin- 
der is cut along the vertical line passing through the center of circle 0 
and “unrolled” (Figure 4.11). The two parastichies and the circumfer- 
ence of the cylinder passing through point 0 form a triangle as shown 
in Figure 4.12. The perpendicular to the base from point mn divides 
this triangle into two right triangles. If d denotes the diameter of the 
circles, then 


(ndm)? + (mnh)? = (nd)? 


and 
(m6,)? + (mnh)? = (md)?. 


The above system of equations can be solved with respect to h and d: 


h = 4/(52, — 63)/(n? —m?) (4.7) 
d = 4/(n262, — m262)/(n? — m?) (4.8) 





or, taking into consideration equation (4.5), 





t= 21 (nbn + m6,,)/(n? — m?). (4.9) 
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Figure 4.13: Patterns of tangent circles drawn on the surface of a cylinder 
as a function of circle diameter 


The problem is to determine values of 6,, and 6,. They are not 
simply functions of parameters m and n. Figure 4.13 shows that, for a 
given m and n, the values of 6,, and 6, can be chosen from a certain 
range, yielding parastichies of different steepness. In order to determine 
this range, observe that at its limits the phyllotactic pattern changes; 
one previously evident parastichy disappears and another is formed. 
Thus, at the range limit, three evident parastichies coexist. It follows 
from Figure 4.13 that at one end of the range the third parastichy has 
order |m—n|, and at the other end it has order (m+n). Three coexisting 
parastichies imply that each circle is tangent to six other circles. In 
other words, all circles lie in the vertices of a regular hexagonal grid, 
as seen in Figure 4.13a and c. Consequently, the angle G+ y at vertex 
mn (Figure 4.12) is equal to 27/3. Expressing the base of the triangle 
in terms of its two sides and their included angle results in 


(27)? = (nd)? + (md)? — 2(nd)(md) cos(27/3) 
or, after simplification, 
d = 2n/Vm? + mn + n?. (4.10) 
Equations (4.9) and (4.10) yield 


NOm + Mb = 2n(n? — m?)/(m? +mn+ n’*). (4.11) 


Pattern 
construction 


Triple-contact 
patterns 
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Solving the system of equations (4.5) and (4.11) with respect to 0. 
and 6, produces 


bm = 1(m + 2n)/(m? + mn +n?) (4.12) 


and 
On = 1(2m +n)/(m? +mn+n?). (4.13) 


Given the values of 6,, and 6,, the divergence angle a can be found 
from either equation (4.3) or (4.4), assuming that the encyclic numbers 
A,» or A,, are known. It follows from the definition that these numbers 
are the smallest positive integers satisfying equation (4.6). A system- 
atic method for solving this equation, based on the theory of continuous 
fractions, is presented by van Iterson [75]. Erickson [36] points out that 
in practice the solution can often be found by guessing. Another pos- 
sibility is to look for the smallest pair of numbers (A,,, A,,) satisfying 
(4.6) using a simple computer program. 

In conclusion, a phyllotactic pattern characterized by a pair of num- 
bers (m,n) can be constructed as follows: 


1. Find A,, and A,, from equation (4.6). 


2. Find the range of admissible values of the angular displacements 
dm and 6,. The limits can be obtained from equations (4.12) and 
(4.13) using the values of m and n for one limit, and the pair 
(min{m,n},|m—n|) for the other. 


3. For a chosen pair of admissible displacement values 6,, and dn, 
calculate the divergence angle a from equation (4.3) or (4.4) and 
the vertical displacement h from equation (4.8). 


4. Find the diameter d of the circles from equation (4.8). 


The diameter d does not enter directly in any formula used for 
image synthesis, but serves as an estimate of the size of surfaces to 
be incorporated in the model. This algorithm was applied to produce 
Table 4.1 showing parameter values for which three parastichies coexist. 
Given a pattern with two parastichies, this table provides the limits of 
the divergence angle a. For example, a (5,8) pattern can be formed for 
values of a ranging from 135.918365° to 138.139542°, which correspond 
to the patterns (3,5,8) and (5,8,13), respectively. 

Further information relating the divergence angle a@ to the verti- 
cal displacement h for various phyllotactic patterns is shown in Fig- 
ure 4.14. The arcs represent parameters of patterns with two paras- 
tichies (m,n). The branching points represent parameters of patterns 
with three parastichies (m,n,m +n). 
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m,n,m+n | a (degrees) h d 





180.000000 | 1.81380 - 
128.571426 | 0.777343 | 2.374821 
96.923073 | 0.418569 | 1.742642 
142.105270 | 0.286389 | 1.441462 
77.142860 | 0.259114 | 1.371104 
102.162163 | 0.147065 | 1.032949 
135.918365 | 0.111049 | 0.897598 
152.307693 | 0.139523 | 1.006115 
63.870968 | 0.175529 | 1.128493 
79.672134 | 0.089203 | 0.804479 
98.709671 | 0.058510 | 0.651536 
107.088600 | 0.068878 | 0.706914 
131.752579 | 0.056097 | 0.637961 
138.139542 | 0.042181 | 0.553204 
150.275223 | 0.049921 | 0.601820 
158.507462 | 0.081215 | 0.767613 
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Table 4.1: Cylinder formula values for triple-contact patterns 











Figure 4.14: The vertical displacement h as a function of the divergence 
angle a for various phyllotactic patterns 


Pineapple 


Spruce cones 


Sedge 


116 Chapter 4. Phyllotaxis 





Figure 4.15: Pineapples 


Models of fruits synthesized using the cylindrical model are shown in 
Figures 4.15 and 4.16. The pineapple (Figure 4.15) is an example 
of a pattern where a given scale has six neighbors, which belong to 
5-, 8- and 13-parastichies. The corresponding divergence angle a is 
equal to 138.139542°. The spruce cones (Figure 4.16) were generated 
using the values m = 5, n = 8 and a = 187.5° (the divergence angle 
a for a (5,8)-parastichy pattern belongs to the interval 135.918365° to 
138.139542°). From these values, h and d were calculated as a function 
of the radius of the cylinder. The effect of closing the bottom and 
top of the pineapple and spruce cones was achieved by decreasing the 
diameter of the cylinder and the size of the scales. 

A variant of the model of phyllotaxis on a cylinder can be used to 
model organs that are conical rather than cylindrical in shape. For 
example, Figure 4.17 shows a model of the sedge Carex laevigata. L- 
system 4.1 generates the male spike. Production p; specifies the basic 
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Figure 4.17: Carex laevigata: the male spike, the entire shoot, the female 
spike 
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#define IRATE 1.025 /* internode growth rate */ 
#define SRATE 1.02 /* spikelet growth rate */ 


#include M /* spikelet shape specification */ 
wi: A 
py: UA 2k — [&(5)£(1)0M(1) ] F(0.2)/(137.5)A 


p2 : MCs): s<3 — M(s*SRATE) 
p3 : £(s): s<3 — £(s*SRATE) 
pa: &(s): s<15— &(s*SRATE) 
ps : FCi): i<t — F(i* IRATE) 


L-system 4.1: The male spike of Carex laevigata 


layout of the spikelets, similar to that given by the L-system in Fig- 
ure 4.11. The top portion of the spike has a conical shape due to the 
growth of spikelets for some time after their creation by the apex. Ac- 
cording to production p2, a spikelet grows by factor SRATE in each 
derivation step, until it reaches the threshold size of 3. In an analogous 
way, productions p3, pa and ps capture the distance increase between 
spikelets and the spike axis, the increase of the branching angle, and 
the elongation of internodes. 

The models of phyllotaxis deal with the arrangement of organs in 
space. For the purpose of mathematical analysis their shape is reduced 
to a simple geometric figure, usually a circle. However, in realistic 
images the exact shape of the organs must be captured as well. Several 
techniques suitable for this purpose are outlined in the next chapter. 


Chapter 5 


Models of plant organs 


Many concepts presented in the previous chapters were illustrated using 
realistic images, but the modeling techniques for leaves and petals have 
not been described yet. Three approaches are discussed below. 


5.1 Predefined surfaces 


The standard computer graphics method for defining arbitrary surfaces 
makes use of bicubic patches [9, 10, 40]. A patch is defined by three 
polynomials of third degree with respect to parameters s and t. The 
following equation defines the x coordinate of a point on the patch: 


x(s,t) = ay, 8°t? + aygs*t? + a438%t + ayas? 
Th anise + asset" + dio387t + dina $* 
+ ag st? + aggst? + aggst + agas 


sk Giayd? + Gag? = dagt aa 





Analogous equations define y(s,t) and z(s,t). All coefficients are de- 
termined by interactively designing the desired shape on the screen 
of a graphics workstation. Complex surfaces are composed of several 
patches. 

The surfaces are incorporated into a plant model in a manner sim- 
ilar to subfigures (Section 1.4.2). The L-system alphabet is extended 
to include symbols representing different surfaces. When the turtle 
encounters such a symbol preceded by a tilde (~), the corresponding 
surface is drawn. 

The exact position and orientation of surface S representing an ap- 
pendage is determined using a contact point Ps, the heading vector Hs 
and the up vector Vs as a reference (Figure 5.1). The surface is trans- 
lated in such a way that its contact point matches the current position 


Bicubic patches 


Turtle 
interpretation 


Examples 


Contour 
tracing 


Fern 
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Pe 


H 


Figure 5.1: Specification of an appendage 


of the turtle, and is rotated to align its heading and up vectors with the 
corresponding vectors of the turtle. If a surface represents an internal 
part of the structure such as an internode, a distinction between the 
entry and exit contact points is made. 

The majority of organs presented in this book have been modeled 
this way. The petals of sunflowers, zinnias, water-lilies and roses shown 
in Chapter 4 provide good examples. Figure 5.2 illustrates an additional 
improvement in the appearance of organs, made possible by the appli- 
cation of textures to the surfaces of leaves, flowers and vine branches. 


5.2 Developmental surface models 


Predefined surfaces do not “grow.” String symbols can be applied to 
control such features as the overall color and size of a surface, but 
the underlying shape remains the same. In order to simulate plant 
development fully, it is necessary to provide a mechanism for changing 
the shape as well as the size of surfaces in time. One approach is to 
trace surface boundaries using the turtle and fill the resulting polygons. 
A sample L-system is given below: 


w: DL 
pp: LDo{-FX+X—-FX—-|-FX+X+4+FxX} 
po: X > PX 


Production p, defines leaf L as a closed planar polygon. The braces 
{ and } indicate that this polygon should be filled. Production p2 
increases the lengths of its edges linearly. The model of a fern shown 
in Figure 5.3 incorporates leaves generated using this method, with 
the angle increment equal to 20°. Note the phase effect due to the 
“growth” of polygons in time. A similar approach was taken to generate 
the leaves, flowers and fruits of Capsella bursa-pastoris (Figure 3.5 on 
page 74). 
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Figure 5.2: Maraldi figure by Greene [54] 





Figure 5.3: The fern 
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4 12 
a b 
9,10 
3 5 
6,7 
2 6 
1 7 1,2,3,4,5,8,11 \ 
Starting point Starting point 


a{I+H++G . |[++GG . |[+GGG , ||GGGGG . ]|-GGG , ]|-GG . | —G .} 
b[+[+[+1+A{.].C.}{.].GC.}{.].GGC.}{.]. GGGC.} 
1 2 3 4 5 6 7 8 9 10 11 12 


Figure 5.4: Surface specification using a tree structure as a framework 


In practice, the tracing of polygon boundaries produces acceptable ef- 
fects only in the case of small, flat surfaces. In other cases it is more 
convenient to use a tree structure as a framework. Polygon vertices are 
specified by a sequence of turtle positions marked by the dot symbol 
(.). An example is given in Figure 5.4a. The letter G has been used 
instead of F to indicate that the segments enclosed between the braces 
should not be interpreted as the edges of the constructed polygon. The 
numbers correspond to the order of vertex specification by the turtle. 

Figure 5.5 shows the development of a cordate leaf modeled using 
this approach. The axiom contains symbols A and B, which initiate the 
left-hand and right-hand sides of the blade. Each of the productions p; 
and pz creates a sequence of axes starting at the leaf base and gradually 
diverging from the midrib. Production p3 increases the lengths of the 
axes. The axes close to the midrib are the longest since they were 
created first. Thus, the shape of this leaf is yet another manifestation 
of the phase effect. The leaf blade is defined as a union of triangles 
rather than a single polygon. Such triangulation is advantageous if 
the blade bends, for example due to tropism (Chapter 2). Figure 5.4b 
provides an additional illustration of the model by magnifying the left 
side of the leaf after four derivation steps. 
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Do : B— [-B 
p3 : C — GC 


Figure 5.5: Developmental sequence of a cordate leaf generated using an 
L-system 


The described method makes it possible to define a variety of leaves. 
Their shapes depend strongly on the growth rates of segments. For 
example, a family of simple leaves and the corresponding parametric 
L-system are shown in Figure 5.6. 

According to production p;, in each derivation step apex A(t) ex- 
tends the main leaf axis by segment G(LA, RA) and creates a pair of 
lateral apices B(t). New lateral segments are added by production po. 
Parameter t, assigned to apices B by production p;, plays the role of 
“erowth potential” of the branches. It is decremented in each derivation 
step by a constant PD, and stops production of new lateral segments 
upon reaching 0. Segment elongation is captured by production p3. 

For the purpose of analysis, it is convenient to divide a leaf blade 
into two areas. In the basal area, the number of lateral segments is 
determined by the initial value of growth potential ¢ and constant PD. 
Since the initial value of t assigned to apices B increases towards the 
leaf apex, the consecutive branches contain more and more segments. 
On the other hand, branches in the apical area exist for too short a 
time to reach their limit length. Thus, while traversing the leaf from 
the base towards the apex, the actual number of segments in a branch 
first increases, then decreases. As a result of these opposite tendencies, 
the leaf reaches its maximum width near the central part of the blade. 


Simple leaves 
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n=20, 6=60° 


#define LA 5 
#define RA 1 
#define LB 1 
#define RB 1 
#define PD 1 


w : {.A(0)} 
Pi: A(t) 

p2 : Bt) 

p3 : G(s,r) 





/* 
/* 
/* 
/* 
/* 


* 
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initial length - main segment */ 
growth rate - main segment */ 
initial length - lateral segment */ 
growth rate - lateral segment */ 
growth potential decrement */ 


— G(LA,RA) [-B(t).] [A(t+1)] [+B(t) .] 


t>O0 — G(LB,RB)B(t-PD) 


* 


— G(s*r,r) 


Figure 5.6: A family of simple leaves generated using a parametric L-system 
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LA|;RA};} LB} RB} PD 
5} 1.0} 1.0} 1.00 | 0.00 
5} 1.0} 1.0} 1.00 | 1.00 
5} 1.0} 0.6} 1.06 | 0.25 
5 | 1.2} 10.0} 1.00 | 0.50 
5} 1.2) 4.0) 1.10 | 0.25 
5} 1.1] 1.0} 1.20 | 1.00 





Table 5.1: Values of constants used to generate simple leaves 
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Figure 5.7: A rose in a vase 


Table 5.1 lists the values of constants corresponding to particular shapes. 


For a given derivation length, the exact position of the branch with the 
largest number of segments is determined by PD. If PD is equal to 0, 
all lateral branches have an unlimited growth potential, and the basal 
part of the leaf does not exist (Figure 5.6a). If PD equals 1, the basal 
and apical parts contain equal numbers of lateral branches (Figures 5.6 
b and f). Finer details of the leaf shape are determined by the growth 
rates. If the main axis segments and the lateral segments have the same 
growth rates (RA = RB), the edges of the apical part of the leaf are 
straight (Figures 5.6 a and b). If RA is less than RB, the segments 
along the main axis elongate at a slower rate than the lateral segments, 
resulting in a concave shape of the apical part (Figures 5.6 ¢ and f). In 
the opposite case, with RA greater than RB, the apical part is convex 
(Figures 5.6 d and e). The curvature of the basal edges can be analyzed 
in a similar way. 


Shape control 
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n=25, d=60° 

#define LA 5 /* initial length - main segment */ 

#define RA 1.15 /* growth rate - main segment */ 

#define LB 1.3 /* initial length - lateral segment */ 

#define RB 1.25 /* growth rate - lateral segment */ 

#define LC 3 /* initial length - marginal notch */ 

#define RC 1.19 /* growth rate - marginal notch */ 

w : [{ACO,0).}] [{ACO,1).}] 

pi: ACt,d) : d=O0 -—  .G(LA,RA).[+B(t)G(LC,RC,t).}] 
[+B(t){.]A(t+1,d) 

p2 : A(t,d) : d=l1 — .G(LA,RA).[-B(t)G(LC,RC,t).}] 
[-B(t){.]A(t+1,d) 

p3 : B(t) : t>0 — G(LB,RB)B(t-1) 

pa : G(s,r) : Ox — G(s*r,r) 


ps : G(s,r,t) : tol — G(s*r,r,t-1) 


Figure 5.8: A rose leaf 


Figure 5.7 shows a long-stemmed rose with the leaves modeled accord- 
ing to Figure 5.8. The L-system combines the concepts explored in 
Figures 5.5 and 5.6. The axiom contains modules A(0,0) and A(0, 1), 
which initiate the left-hand and right-hand side of the leaf. The de- 
velopment of the left side will be examined in detail. According to 
production p;, in each derivation step apex A(t,0) extends the midrib 
by internodes G(LA, RA) and creates two colinear apices B(t) pointing 
to the left. Further extension of the lateral axes is specified by produc- 
tion p3. The leaf blade is constructed as a sequence of trapezoids, with 
two vertices lying on the midrib and the other two vertices placed at 
the endpoints of a pair of lateral axes formed in consecutive derivation 
steps. The module G(LC, RC,t) introduces an offset responsible for 
the formation of notches at the leaf margin. Production p4 describes 
the elongation of internodes responsible for overall leaf shape, while 
production ps controls the size of the notches. The development of the 
right side of the blade proceeds in a similar manner, with production po 
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Figure 5.9: Surface specification using stacked polygons 


recreating the midrib and creating lateral apices pointing to the right. 
The bending of the midrib to the right is a result of tropism. 

In the examples discussed so far, the turtle specifies the vertices of 
one polygon, then moves on to the next. Further flexibility in surface 
definition can be achieved by interleaving vertex specifications for dif- 
ferent polygons. The turtle interpretation of the braces is redefined in 
the following way. A string containing nested braces is evaluated us- 
ing two data structures, an array of vertices representing the current 
polygon and a polygon stack. At the beginning of string interpretation, 
both structures are empty. The symbols {, } and . are then interpreted 
as follows: 


{ Start a new polygon by pushing the current polygon on 
the polygon stack and creating an empty current polygon. 


Append the new vertex to the current polygon. 


} Draw the current polygon using the specified vertices, 
then pop a polygon from the stack and make it the current 


polygon. 


An example of string interpretation involving nested braces is given in 
Figure 5.9. 

The above technique was applied to construct the flowers of the 
lily-of-the-valley shown in Figure 3.4 (page 72), and magnified in Fig- 
ure 5.10. A flower is represented by a polygon mesh consisting of five 
sequences of trapezoids spread between pairs of curved lines that em- 
anate radially from the flower base. A single sequence is generated by 
the following L-system: 


w : [X(36)A]/(72)[X(36)B] 


pi: A * — [&GA{.]. 
po: Bo: x > B&.G.} 
p3: X(a): * — X(a+4.5) 


Nested 
polygons 
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valley 
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branching 
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Figure 5.10: Structure of a lily-of-the-valley flower 


Productions p; and ps create two adjacent framework lines and mark 
polygon vertices consistently with Figure 5.9. Production ps controls 


the angle at which the framework lines leave the stalk at the flower 
base. 


5.3. Models of compound leaves 


So far, the discussion of organ models has focused on the definition 
of surfaces. However, in the case of highly self-similar structures, the 
individual surfaces become inconspicuous, and the expression of the 
geometric relationships between younger and older parts of the struc- 
ture becomes the key issue. For example, Figure 5.11 shows compound 
leaves often found in the family Umbelliferae. According to produc- 
tion pz, the apex A(0) creates two segments F'(1) and a pair of lateral 
apices A(D) in each derivation step. Production p; delays the develop- 
ment of the daughter branches by D steps with respect to the mother 
branch. This pattern is repeated recursively in branches of higher or- 
der. Production p3 gradually elongates the internodes, and in this way 
establishes proportions between parts of a leaf. The values of the con- 








Figure | D| R_ | Derivation length 
a 0 | 2.00 10 
b 1 | 1.50 16 
c 2 | 1.36 21 
d 4 | 1.23 30 
e 7) 1.17 45 














Table 5.2: Values of constants used to generate compound leaves 


5.3. Models of compound leaves 


A eas 13 
PAF VY 
5 \ Be 





4 i 
Oe) v% i 
oy 4 


#define D i 


#define R 1.5 


W 
Pi 


p2 : 
P3 


A(O) 
A(d) 
A(d) 
F(a) 


/* 
/* 


*aQad 
ll Vv 
oo 


Figure 5.11 
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#define D i 
#define R 1.36 
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/* apical delay */ 
/* internode elongation rate */ 


A(d-1) 
F(1) [+A(D)] F(1)B(0) 
B(d-1) 
F(1) [-B(D) ] F(1) A(O) 
F (a*R) 
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Figure 5.12: Compound leaves with alternating branching patterns 








Figure | D Derivation length 
a 1 | 1.36 20 
b 4 | 1.18 34 
c #| del3 46 














Table 5.3: Values of the constants used to generate compound leaves with 
alternating branches 
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stants used to generate the structures shown in Figure 5.11 are listed 
in Table 5.2. The model is sensitive to growth rate values — a change 
of 0.01 visibly alters proportions. Leaves of the wild carrot, shown in 
Figure 3.23 (page 96), correspond closely to case Figure 5.11b. 

Another type of compound leaf, with alternating branches, is ex- 
amined in Figure 5.12. The values of the constants used in the L- 
system are specified in Table 5.3. Further examples can be found in 
the next chapter, devoted to the animation of developmental processes, 
and Chapter 8, which ties developmental models with fractals. 


Alternating 
branching 


Chapter 6 


Animation of plant 
development 


The sequences of images used in Chapters 3 and 5 to illustrate the de- 
velopment of inflorescences and compound leaves suggest the possibil- 
ity of using computer animation to visualize plant development. From 
a practical perspective, computer animation offers several advantages 
over traditional time-lapse photography. 


Photography is sensitive to imperfections in the underlying ex- 
periment. A disease or even a temporary wilting of a plant may 
spoil months of filming. 


In nature, developmental processes are often masked by other 
phenomena. For example, the growth of leaves can be difficult 
to capture because of large changes in leaf positions during the 
day. Similarly, positions of tree branches may be affected by 
wind. Computer animation makes it possible to abstract from 
these distracting effects. 


Animation can be used when time-lapse photography is imprac- 
tical because of the long development time of plants (e.g. years 
in the case of trees). 


In time-lapse photography, the initial position of the camera and 
the light conditions must be established a priori, without know- 
ing the final shape of the plant. In computer animation all de- 
velopmental stages of the modeled plant are known in advance, 
allowing for optimal selection of the model view. 


Animation can be applied to visualize developmental mechanisms 
that cannot be observed directly in real plants, such as the prop- 
agation of hormones and nutrients. 


Motivation 


Discrete 
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L-systems 
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e Animation offers an unprecedented means for visualizing the de- 
velopment of extinct plants on the basis of paleobotanical data. 


The original formalism of L-systems provides a model of develop- 
ment that is discrete both in time and space. Discretization in time 
implies that the model states are known only at specific time inter- 
vals. Discretization in space means that the range of model states is 
limited to a finite number. Parametric L-systems remove the latter ef- 
fect by assigning continuous attributes to model components. However, 
the model states are still known only in discrete time intervals. This 
presents a problem in animation, where a smooth progression of the 
developing forms is desirable. 

This last statement requires further clarification. The very nature 
of animation is to produce the impression of continuous motion by dis- 
playing a sequence of still frames, captured at fixed time intervals. Why 
is a continuous model of plant development needed if it will be used to 
generate a fixed sequence of images in the final account? Wouldn’t it 
be enough to retain the standard definition of L-systems and assume 
time slices fine enough to produce the desired progression of forms? 
This approach, while feasible and useful, has three major drawbacks. 


e A model can be constructed assuming longer or shorter time in- 
tervals, but once the choice has been made, the time step is a part 
of the model and cannot be changed easily. From the viewpoint 
of computer animation it is preferable to control the time step by 
a single parameter, decoupled from the underlying L-system. 


e The continuity criteria responsible for the smooth progression of 
shapes during animation can be specified more easily in the con- 
tinuous time domain. 


e It would be conceptually elegant to separate model development, 
defined in continuous time, from model observation, taking place 
in discrete intervals. 


A developmental process is viewed as consisting of periods of contin- 
uous module expansion delimited by instantaneous module divisions. 
Special conditions are imposed to preserve the shape and growth rates 
of the organism during these qualitative changes. An analogy can be 
drawn to the theory of morphogenesis advanced by Thom [142], who 
viewed shape formation as a piecewise continuous process with singu- 
larities called catastrophes. 

Formally, development taking place in continuous time is captured 
using the formalism of timed DOL-systems. The key difference between 
these L-systems and the types of L-systems considered so far lies in 
the definition of the derivation function. In “ordinary” L-systems, the 
derivation length is expressed as the number of derivation steps. In 
timed DOL-systems, the derivation length is associated with the total 
elapsed time since the beginning of the observation. 
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6.1 Timed DOL-systems 


Let V be an alphabet and R the set of positive real numbers (including 
0). The pair (a,7) € V x R is referred to as the timed letter a, and 
the number 7 is called the age of a. A sequence of timed letters, x = 
(41,71) ---(@n,T) € (V x R)*, is called a timed word over alphabet V. 

A timed DOL-system (tDOL-system) is a triplet G = (V,w,P), 
where 


e V is the alphabet of the L-system, 


e w €(V x R)* is a nonempty timed word over V, called the initial 
word, and 


e PC(V x R) x (V x R)* is a finite set of productions. 


Instead of writing ((a, (3), (b1,01)...(bn,@n)) € P, the notation 
(a, 3) — (b1,01)...(bn,Qn) is used. The parameter ( is referred to 
as the terminal age of the letter a, and each parameter a; is the initial 
age assigned to the letter b; by production P. The following assump- 
tions are made: 


Cl. For each letter a € V there exists exactly one value 3 € R such 
that (a, 3) is the predecessor of a production in P. 


C2. If (a,@) is a production predecessor and (a,q) is a timed let- 
ter that occurs in the successor of some production in P, then 
B>a. 


According to these conditions, each letter has a uniquely defined 
terminal age. Furthermore, an initial age assigned to a letter by a 
production must be smaller than its terminal age, i.e., its lifetime (G — 
a) must be positive. 

Let (a, 3) — (b1,01)...(bn,Q@,) be a production in P. A function 
D: ((V x R)* x R) > (V x R)* is called a derivation function if it has 
the following properties: 


Al. D(((a1, Pies - (Ga, Tn))s t) = D((a, a); t) auens D((an, Tn); t) 

A2. D((a,7),t) = (a,7 +1), if r+t< p 

A3. D((a,7),t) = D((b1, 01)... (bn, An),t -(G—T)), if T+t> 

A derivation in a timed DOL-system is defined in terms of two types 


of time variables. Global time ¢t is common to the entire word under 
consideration, while local age values 7; are specific to each letter aj. 


Definition 


Derivation 
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Axiom A1 identifies t as the variable that synchronizes the entire de- 
velopment, and specifies that the lineages of all letters can be considered 
independently from each other (thus, no interaction between letters is 
assumed). With the progress of time t, each letter “grows older” until 
its terminal age is reached (axiom A2). At this moment subdivision 
occurs and new letters are produced with initial age values specified by 
the corresponding production (axiom A3). Condition C1 guarantees 
that the subdivision time is defined unambiguously, hence the devel- 
opment proceeds in a deterministic fashion. Condition C2 guarantees 
that, for any value of time t, the recursive references specified by axiom 
A3 will eventually end. 

The above concepts are examined by formulating a timed DOL- 
system that simulates the development of a vegetative part of the An- 
abaena catenula filament. Given the discrete model expressed by equa- 
tion (1.1) on page 5, the corresponding tDOL-system is as follows: 


Ww: (a,,0) 

Pi: (a,, 1) = (ai, 0)(b,, 0) 

P2: (az, 1) a (bi, 0)(a,, 0) (6.1) 
P3 : (Os; 1) a (a,, 0) 

Pa: (by, 1) ae (ai, 0) 


In accordance with the discrete model, it is assumed that all cells have 
the same lifetime, equal to one time unit. The derivation tree is shown 
in Figure 6.1. The nodes of the tree indicate production applications 
specified by axiom A3, and the triangular “arcs” represent the continu- 
ous aging processes described by axiom A2. The vertical scale indicates 
global time. For example, at time t = 2.75 the filament consists of three 
cells, b;, a, and a,, whose current age is equal to 0.75. 

According to the definition of time intervals corresponding to axioms 
A2 and A38, a production is applied after the age 7 + t exceeds the 
terminal age. Consequently, at division time the “old” cells still exist. 
For example, at time t = 2.0 the filament consists of two cells, a; and 
b,, both of age 7 = 1. 

The above L-system can be simplified by considering cells of type b 
as young forms of the cells of type a. This is suggested by Figure 6.1 
where cells b simply precede cells a in time. The simplified L-system 
has two productions: 


es Ces 2) Sap Cae 
os ar = ea ia (6.2) 


The corresponding derivation tree starting from cell (a,,1) is shown in 
Figure 6.2. Note the similarity to the tree from the previous example. 

Whether a natural developmental process or its mathematical model 
is considered, the choice of observation times and the act of observation 
should not affect the process itself. In other terms, each derived word 
should depend only on the total elapsed time t, and not on the partition 
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Time 
0.0 


2.0 





2.75 
3.0 


4.0 





5.0 
Figure 6.1: Derivation tree representing the continuous development of An- 


abaena catenula described by the L-system is equation (6.1). Sections of the 
triangles represent cell ages. 


Time 
0.0 


2.0 





215 
3.0 


4.0 





5.0 


Figure 6.2: Derivation tree representing the continuous development of An- 
abaena catenula corresponding to the rules specified in equation (6.2) 


138 Chapter 6. Animation of plant development 


of ¢ into intervals. The following theorem shows that timed DOL- 
systems satisfy this postulate. 

Theorem. Let G = (V,w, P) be atDOL-system, and x € (V x R)* be 
a timed word over V. For any values of t,,t, > 0, the following holds: 


D(D(z, ta), i) = D(z, ta + te) 


Proof. Let us first consider the special case where the word x consists 
of a single timed letter, (a9, 7), and all productions in the set P take 
single letters into single letters. According to condition C1, there exists 
a unique sequence of productions from P such that: 


Rae Oe) Gehiga tac = Oyo ci 


Let (a%,7%) be the result of the derivation of duration t, that starts 
from (ao,7). According to axioms A2 and A3, and assuming that 
ta > 29 — 7, this derivation can be represented in the form 


D ((ao, To); ta) 
= D((a1,01), ta — (Bo —70)) 
= D((a2, a2), ta — (G0 — 70) — (G1 — O1)) 


= D((ay, Ox), ta — (G0 — 7) — (G1 — G1) —- +. = (Be_-1 — On-1)) 


= (ax, Tr), 
where 
yee 
Th = On + [ta — (Bo — 70) — 95 (8; — a4)]. 
1=1 


Since the sequence of recursive calls can be terminated only by an ap- 
plication of axiom A2, the index k and the age 7; satisfy the inequality 


Op TE SS De 


Due to condition C2, such an index k always exists and is unique. 

Let us now consider a derivation of duration ty > 9, — T, that 
starts from (az,7%). Following the same reasoning, the result can be 
represented as (@m,T7m), Where 


m—1 


Tm = Om + [ty — (Bk — Th) — » (2; — a%)] 


i=k+1 
and 
Oa = Ta = One: 
By substituting the value of 7; into the formula for 7,,, we obtain after 
simple transformations 


—1 


Tm = Om + [(ta + te) — (Go — 70) — D_, (fi — 04)]. 


w=1 


3 
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Thus, the timed letter (Qj, 7m) also results from the derivation of du- 
ration t, +t, starting with (ag, 79): 


(Qm, Tm) = D(D((40, 70), ta); to) = D((G0, 70), ta + ty). 
So far, we have considered only the case 
ta > Bo-—T, te > Be—- Tr. 
Two other cases, namely, 
O<ta<Po-—7%, ty > Be-Te 
and 
ta > Po-—T, Oto < Be—T 
can be considered in a similar way. The remaining case, 
OStS fo—t Ofte Ss Be — te; 


is a straightforward consequence of condition C2. This completes the 
proof of the special case. In general, a derivation that starts from a 
word (@1,71)---(@n,;7) can be considered as n separate derivations, 
each starting from a single letter. This observation applies not only to 
the initial word specified at time t = 0, but also to any intermediate 
word generated during the derivation. Consequently, any path in the 
derivation tree can be considered as a sequence of mappings that takes 
single letters into single letters. Application of the previous reasoning 
separately to every path concludes the proof. O 





6.2 Selection of growth functions 


Timed L-systems capture qualitative changes in model topology corre- 
sponding to cell (or, in general, module) divisions, and return the age 
of each module as a function of the global time t. In order to complete 
model definition, it is also necessary to specify the shape of each module 
as a function of its age. Potentially, such growth functions can be es- 
timated experimentally by observing real organisms [72, 73]. However, 
if detailed data is not available, growth functions can be selected from 
an appropriate class by choosing parameters so that the animation is 
smooth. This approach can be viewed as more than an ad hoc tech- 
nique for constructing acceptable animated sequences. In fact, Thom 
presents it as a general methodology for studying morphogenesis [142, 
page 4]: 


The essence of the method to be described here consists 
in supposing a priori the existence of a differential model 
underlying the process to be studied and, without knowing 
explicitly what the model is, deducing from the single as- 
sumption of its existence conclusions relating to the nature 
of the singularities of the process. 


Continuity 
requirements 


Linear growth 
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A technique for computing parameters of growth functions in the 
case of nonbranching filaments and simple branching structures is given 
below. 


6.2.1 Development of nonbranching filaments 


In a simple case of geometric interpretation of timed L-systems, symbols 
represent cells that elongate during their lifetime and divide upon reach- 
ing terminal age. Model geometry does not change suddenly, which 
means that 


e the length of each cell is a continuous function of time, and 


e the length of a cell before subdivision is equal to the sum of the 
lengths of the daughter cells. 


The above continuity requirements are formalized in the context of a 
tDOL-system G = (V,w, P) as follows: 


R1. Let [@min, @] describe the life span of a timed letter a € V. The 
age Amin is the minimum of the initial age values assigned to a 
by the axiom w or by some production in P. The terminal age 
G is determined by the predecessor of the production acting on 
symbol a. The growth function g(a,7), which specifies the length 
of cell a as a function of age 7, must be a continuous function of 
parameter 7 in the domain [Qmin, (I. 


R2. For each production (a, 3) — (61,01)... (0n,Qn) in P the follow- 
ing equality holds: 


n 


g(a, 3) = S> 9(bi, a) (6.3) 


w=1 


In practice, requirement R1 is used to select the class of growth func- 
tions under consideration, and the equations resulting from requirement 
R2 are used to determine the parameters in function definitions. 

For example, in the case of the timed L-system specified in equa- 
tion (6.2), requirement R2 takes the form 


g(ar,2) = g(a,1)+ g(a,,0) 
g(a1,2) = g(a,0) + g(ar, 1). (6.4) 


Let us assume that the growth functions are linear functions of time: 


g(ai,7T) = AT+B, 


g(ar,T) = A,7+ B, (6.5) 
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Figure 6.3: Diagrammatic representation of the development of Anabaena 
catenula, with (a) linear and (b) exponential growth of cells 


By substituting equations (6.5) into (6.4), we obtain 





2A, + B, = (1A, Bi) (0A, B,) or 2A, = A, + B, 
2A, + B, = (0A; + Bi) (1A, Ba) or 2A, = A, + B,. 











The desired continuity of development is provided by all solutions of 
this system. They can be expressed in terms of coefficient c, which 
relates the growth rate of cells a; to that of cells a;,: 


Aj = cA, 
B, = A,(2 = C) 
B, = A,(2c—1) 


Figure 6.3a illustrates the developmental process considered for c = 
1. The diagram is obtained by plotting the cells in the filament as 
horizontal line segments on the screen. Colors indicate cell type and 
age. The observation time t ranges from 1 (at the top) to 7 (at the 
bottom), with increment At = 0.009. 

The slopes of the side edges of the diagram represent growth rates 
of the entire structure. Notice that they remain constant in the peri- 
ods between cell divisions, then change. This effect is disconcerting in 
animation, since the rate of organism growth suddenly increases with 
each cell division. In order to prevent this, it is necessary to extend re- 
quirements R1 and R2 to a higher order of continuity N. Specifically, 
equation (6.3) takes the form 


=D (b:0,,) tor k= 0, Lye. 5, (6.6) 


Exponential 
growth 
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where g*)(a,7) is the k*" derivative of the growth function g(a, 7) with 
respect to age T. 

In the case of Anabaena, an attempt to achieve first-order continu- 
ity assuming linear growth functions yields an uninteresting solution, 
g(a,T) = g(a,,T) = 0. Thus, more complex growth functions must be 
considered, such as an exponential function that can be used to approx- 
imate the initial phase of sigmoidal growth. Assume that the growth 
function has the form 


g(a,T) = g(ap,T) = AeP*. (6.7) 


The objective is to find values of parameters A and B that satisfy 
equation (6.6) for k = 0,1. By substituting equation (6.7) into (6.6), 
we obtain 

AB*e?? — AB*e® + AB, (6.8) 


which implies that zero-order continuity yields continuity of infinite 
order in this case. Solution of equation (6.8) for any value of k yields 


1475 
9 


B=ln 





~ 0.4812. Go) 


Parameter A is a scaling factor and can be chosen arbitrarily. The 
corresponding diagrammatic representation of development is shown in 
Figure 6.3b. The side edges of the diagram, representing the growth 
rates of the whole structure, are smooth exponential curves. 


6.2.2 Development of branching structures 


The notions of tDOL-system and growth function extend in a straight- 
forward way to L-systems with brackets. For example, the following 
tDOL-system describes the recursive structure of the compound leaves 
analyzed in Section 5.3. 


w : (a,0) 
pr: (a,1) — (s,0)[(b, 0)][(, 0)](a, 0) 
pr: (b,8) — (a,0) 


According to production p;, apex a produces an internode s, two lateral 
apices b and a younger apex a. Production p2 transforms the lateral 
apices b into a after a delay 3. The daughter branches recursively repeat 
the development of the mother branch. 

Let us assume that the leaf development is first-order continuous, 
yielding the following equations for k = 0, 1: 


(30) +9 (a,0) (6.10) 
(6.11) 


g 
g(b,0) = 0 
a") (a,0) (6.12) 


g'(b, B) = 
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Equation (6.10) provides continuity along an existing axis upon the 
formation of a new internode s. Equation (6.11) specifies that a newly 
formed lateral apex b has zero length and zero growth rate. Equa- 
tion (6.12) guarantees smooth transformation of apices b into apices a. 

Let us further assume that apices a and internodes s expand expo- 
nentially, 


Gast) = Aye?s" (6.13) 
g(s,T) = Ase”, (6.14) 


for 7 € [0,1]. The expansion of the lateral apices cannot be described 
by an exponential function, since it would not satisfy equations (6.11). 
Consequently, a polynomial growth function g(b,7) is chosen. Equa- 
tions (6.11) and (6.12) fix the function’s endpoints and the tangents at 
the endpoints. Thus, in general, g(b,7) must be a polynomial of degree 
three or more, such as 


g(b, 7) = Apr® + Ber? + Cyr + Dp. (6.15) 


The system of equations (6.10) through (6.12) is solved using the initial 
size A, and the growth rate coefficient B, as independent variables. 
By substituting (6.13) and (6.14) into equations (6.10) for k = 0,1, we 
obtain 


Ay = Ayfere a1) 
Bea B.. 


Equations (6.11) and (6.15) yield 
Cy, = D, = 0. 
Finally, substitution of (6.13) and (6.15) into (6.12) results in 


Ap = i (Be — 2) 


Ag 
B= Be 3 — bBo). 
Figure 6.4 shows a sequence of images produced by this model using 
values G = A, = 1 and B, = 0.48. The branching angles are equal to 
45°. The observation time t ranges from 6.9 to 7.7, with an increment 
of 0.2. Note the gradual formation of lateral segments. 

In the examples considered above, modules are represented as 
straight lines, with growth functions controlling line lengths. Other 
parameters, such as the branching angle, the diameter of segments and 
the size of predefined surfaces, can be controlled in an analogous way. 
Generally, any developmental model captured by an OL-system with 
turtle interpretation can be converted into a tDOL-system and ani- 
mated. 
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Figure 6.4: Developmental sequence of a branching structure modeled using 
a tDOL-system 





Chapter 7 


Modeling of cellular layers 


String L-systems, the first formalism considered in this book, are suit- 
able for the modeling of nonbranching filaments such as Anabaena 
catenula. The introduction of brackets extends the class of modeled 
structures to axial trees. However, many structures found in botany 
have a more complex topology, which can only be described by graphs 
with cycles. The developmental surface models presented in Section 5.2 
make it possible to specify a limited class of these graphs. This chapter 
describes a more general approach and applies it to simulate the devel- 
opment of single-layered cellular structures such as those found in fern 
gametophytes, animal embryos and plant epidermis. All structures con- 
sidered are of microscopic dimensions and relatively undifferentiated, 
yet the presented methods may bring us closer to the modeling of more 
complex patterns, such as the venation of leaves. 

The modeling method consists of two stages. First, the topology 
of the cell division patterns are expressed using the formalism of map 
L-systems, which allows for the formation of cycles in a structure. At 
this stage the neighborhood relations between cells are established, but 
the cell shapes remain unspecified. Next, cell geometry is modeled 
using a dynamic method that takes into account the osmotic pressure 
inside the cells and the tension of cell walls. The development can be 
animated by considering periods of continuous cell expansion, delimited 
by instantaneous cell divisions. 


7.1 Map L-systems 


From a mathematical perspective, cellular layers can be represented 
using a class of planar graphs with cycles, called maps [148]. Nakamura 
et al. [102] characterize them as follows: 


e A map is a finite set of regions. Each region is surrounded by 
a boundary consisting of a finite, circular sequence of edges that 
meet at vertices. 


Maps as 
models of cell 
layers 


Map rewriting 


mBPMOL- 
systems 


Definition 
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e Each edge has one or two vertices associated with it. The one- 
vertex case occurs when an edge forms a loop. The edges cannot 
cross without forming a vertex and there are no vertices without 
an associated edge. 


e Every edge is a part of the boundary of a region. 


e The set of edges is connected. Specifically, there are no islands 
within regions. 


A map corresponds to a microscopic view of a cellular layer. Regions 
represent cells, and edges represent cell walls perpendicular to the plane 
of view. The internal components of a cell are not considered. 

The process of cell division can be described by map rewriting. This 
notion is an extension of string rewriting as discussed in Section 1.2. 
In general, map-rewriting systems are categorized as sequential or par- 
allel, and can be region-controlled or edge-controlled [87]. Since several 
cells may divide concurrently, a parallel rewriting system is needed. 
The second categorization has to do with the form of rewriting rules, 
which may express cell subdivisions in terms of region labels or edge 
labels. Both approaches are suitable for biological modeling purposes 
[22]. This chapter focuses on the edge-controlled formalism of Binary 
Propagating Map OL-system with markers, or mBPMOL-systems. It 
was proposed by Nakamura, Lindenmayer and Aizawa |102] as a refine- 
ment of the basic concept of map L-systems introduced by Lindenmayer 
and Rozenberg [91]. The name is derived as follows. A map OL-system 
is a parallel rewriting system that operates on maps and does not al- 
low for interaction between regions. In other words, regions are modi- 
fied irrespective of what happens to neighboring regions (a context-free 
mechanism). The system is binary because a region can split into at 
most two daughter regions. It is propagating in the sense that the edges 
cannot be erased, thus regions (cells) cannot fuse or die. The markers 
specify the positions of inserted edges that split the regions. 

The choice of mBPMOL-systems as a modeling tool has two jus- 
tifications. First, they are more powerful than other interactionless 
map-rewriting systems described in the literature [19, 22, 23]. In addi- 
tion, markers have a biological counterpart in preprophase bands of mi- 
crotubules, which coincide with the attachment sites for division walls 
formed during mitosis [55]. It should be noted, however, that double 
wall systems, introduced by J. and H. B. Liick [93], may be relatively 
easier to specify [23]. 

An mBPMOL-system G consists of a finite alphabet of edge labels 
i, a starting map w with labels from ™, and a finite set of edge pro- 
ductions P. In general, the edges are directed, which is indicated by a 
left or right arrow placed above the edge symbol. In some cases, the 
edge direction has no effect on the system operation. Such an edge 
is called neutral and no arrow is placed above the symbol denoting it. 
Each production is of the form A — a, where the directed or neutral 
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a b Cc 
A—*DC[-E]BF A—*FB[+E]CD A—»B[-B]x[+B]B 
A A A 

i eS 





Figure 7.1: Examples of edge productions 


edge A € ®& is called the predecessor, and the string a, composed of 
symbols from » and special symbols |,],+, and —, is called the succes- 
sor. The symbols outside square brackets specify the edge subdivision 
pattern. Arrows can be placed above edge symbols to indicate whether 
the successor edges have directions consistent with, or opposite to, the 
predecessor edge. Pairs of matching brackets | and | delimit markers, 
which specify possible attachment sites for region-dividing walls. The 
markers are viewed as short branches that can be connected to form a 
complete wall. The strings inside brackets consist of two symbols. The 
first symbol is either + or —, indicating whether the marker is placed 
to the left or to the right of the predecessor edge. The second symbol 
is the marker label, with or without an arrow. The left arrow indicates 
that the marker is directed towards the predecessor edge, and the right 
arrow indicates that the marker is directed away from that edge. If no 
arrow is present, the marker is neutral. 

For example, in the production A > po[-z]BF, the directed prede- 
cessor A splits into four edges D, C, B and F’, and produces a marker 
(Figure 7.la). Successor edges D and B have the same direction as A, 
edge C' has the opposite direction, and F’ is neutral. Marker EF is placed 
to the right of A and is directed towards A. Note that this same pro- 
duction could be written as 4 — rB[+E]CD (Figure 7.1b). As an example 
of a production with a neutral predecessor, consider 4—B[-B]x[+B)B- In 
this case the result of production application does not depend on the 
assumed direction of the predecessor edge (Figure 7.1c). 


A derivation step in an mBPMOL-system consists of two phases. 


e Each edge in the map is replaced by successor edges and markers 
using the corresponding edge production in P. 


e Each region is scanned for matching markers. 
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w : ABAB 
pi : A —  BL[-A] [+A]B 
p, : B — A 


Figure 7.2: Example of a map L-system. In the first step, a distinction 
is made between the edge-rewriting phase and the connection of matching 
markers. 


Two markers are considered matching if: 
e they appear in the same region, 
e they have the same label, and 


e one marker is directed away from its incident edge while the other 
is directed towards the edge, or both markers are neutral. 


If a match is found, the markers are joined to create a new edge that 
will split the region. The search for matching markers ends with the 
first match found, even though other markers entering the same region 
may also form a match. From the user’s perspective, the system be- 
haves nondeterministically since it chooses the pair of markers to be 
connected. The unused markers are discarded. 

The operation of mBPMOL-systems is illustrated in the following 
examples. The L-system shown in Figure 7.2 has two productions. 
Production p; creates markers responsible for region division, while 
production p2 introduces a delay needed to subdivide the regions alter- 
nately by horizontal and vertical edges. 

The L-system shown in Figure 7.3 is a modified version of the pre- 
vious one. The only difference is the addition of an edge x, which 
separates the markers in the successor of production p;. This edge 
creates a Z-shaped offset between the inserted edges A. Z-offsets and 
symmetric S-offsets (Figure 7.4) can be observed in many biological 
structures [22, 92]. 

Figure 7.5 illustrates the operation of an mBPMOL-system with 
directed edges. Productions p, and p3 create markers. Production p4 
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p, : A —  B[-A]x[+A]B 
po: B — A 


Figure 7.3: Example of a map L-system 


em 


Figure 7.4: Offsets between four regions that result from the division of two 
regions sharing a common edge: (a) Z-offset, (b) S-offset 
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Figure 7.5: Example of a map L-system with directed edges 
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transforms edge D into C’, so that in each derivation step there is a pair 
of edges A and C’ to which productions p, and p3 apply. Production p2 
indicates that edges B do not undergo further changes.! The resulting 
structure is that of a clockwise spiral. 


7.2 Graphical interpretation of maps 


Maps are topological objects without inherent geometric properties. In 
order to visualize them, some method for assigning geometric interpre- 
tation must be applied. In its description, the biologically-motivated 
terms cell and wall will be used instead of their mathematical counter- 
parts, region and edge. 

Siero, Rozenberg and Lindenmayer [131] proposed a method that, 
in the simplest case, is expressed by the following rules: 


e walls are represented by straight lines, 


e the starting map is represented by a regular polygon, bounded by 
the walls specified in the axiom, 


e when a production subdivides a wall, all successor walls are of 
equal length, and 


e the position of a wall resulting from the union of two matching 
markers is based on the position of these markers. 


This wall subdivision method was used to draw Figures 7.2, 7.3 and 
7.5. However, in a biological context it creates cells with shapes that 
are seldom observed in nature. 

De Does and Lindenmayer [24] proposed a center of gravity method 
that produces more realistic shapes. According to this method each 
interior vertex of the map is placed at the center of gravity of its neigh- 
bors. Such positioning of vertices has a sound biological justification; 
it minimizes hypothetical forces acting along cell walls, thus bringing 
the entire structure to a state of minimum energy. However, if all ver- 
tices were positioned this way, the entire structure would collapse. To 
counteract this effect, the vertices on the map perimeter are displaced 
outward a fixed distance. Unfortunately, this lacks biological justifi- 
cation and introduces sudden shape changes after cell divisions have 
occurred, making it unsuitable for animation purposes. 

Assuming a dynamic point of view, the shape of cells and thus the 
shape of the entire organism results from the action of forces. The 
unbalanced forces due to cell divisions cause the gradual modification 
of cell shapes until an equilibrium is reached. At this point, new cell 
divisions occur, and the search for an equilibrium resumes. 





‘In further L-systems such identity productions are omitted. 
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The dynamic method is based on the following assumptions: 
e the modeled organism forms a single cell layer, 


e the layer is represented as a two-dimensional network of masses 
corresponding to cell vertices, connected by springs that corre- 
spond to cell walls, 


e the springs are always straight and obey Hooke’s law, 


e each cell exerts pressure on its bounding walls directly propor- 
tional to the wall length and inversely proportional to the cell 
area, 


e the pressure on a wall is divided evenly between the wall vertices, 
e the motion of masses is damped, and 


e no other forces (for example, due to friction or gravity) are con- 
sidered. 


The position of each vertex, and thus the shape of the layer, is com- 
puted as follows. As long as an equilibrium is not reached, unbalanced 
forces put masses in motion. The total force #, acting on a vertex X 
is given by the formula 


Fp = S° Fy te Fi, 
wEew 


where fF, are forces contributed by the set W of walls w incident to 
X, and Fy = —bv is a damping force, expressed as the product of a 
damping factor b and vertex velocity v. 

A wall w € W contributes three forces acting on X (Figure 7.6). 
The tension F, acts along the wall, and its magnitude is determined by 
Hooke’s law, . 

F, = —k(l— lp), 


where k is the spring constant, / is the current spring length, and Ip is 
its rest length. The remaining forces, P, and Pz, are due to the pressure 
exerted by the cells on the left and right sides of the wall. Each force 
acts in a direction perpendicular to the wall and is distributed equally 
between the incident vertices. The magnitude of these forces equals 
p-l, where p is the internal cell pressure and / is the wall length. 
The pressure is assumed to be inversely proportional to the cell area: 
p ~ A‘. This assumption is derived from the equation describing 
osmotic pressure, p = SRT, as a function of the concentration of the 
solute S (n moles per volume V), the ideal gas constant R, and the 
absolute temperature T [129]. Assuming that the cell volume V is 
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Figure 7.6: Forces acting on a cell vertex X 


proportional to the area A captured by the two-dimensional model 
under consideration (V = Ah), pressure can be expressed as 


nRT 
Ah — 
Thus p ~ A“, provided that the term nRT/h is constant. A convenient 
formula for calculating the area A is 


A= [90 (i — tis) (yi + Yi) /21, 


i=1 
where (x;, y;) are coordinates of the M vertices of the cell, ry41 = 21, 
and ym4i = Y1- 
The force F, acts on a mass placed at a map vertex. Newton’s 
second law of motion applies, 
ae. a 
m— = Fr, 
di 7 
where Z is the vertex position. If the entire structure has N vertices, a 
system of 2N differential equations is obtained, 





du; = 
mMm— = Fr, (£1,---, EN, 4; 
a 1, (£4 N; Ui) 
dx; 3 
= Vi; 
dt 
where 7 = 1,2,...,N. The task is to find the sequence of positions 
£1,...,€n at given time intervals, assuming that the functions F;, and 
the initial values of all variables 7?,...,£% and vp,..., 0) are known. 


These initial values are determined as follows. 
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e Coordinates of the vertices of the starting map are included in 
the input data for the simulation. 


e Positions of existing vertices are preserved through a derivation 
step. New vertices partition the divided walls into segments of 
equal length. The initial velocities of all vertices are set to zero. 


The system of differential equations with the initial values given 
above represents an initial value problem. It can be solved numerically 
using the forward (explicit) Euler method |41]. To this end, the dif- 
ferential equations are rewritten using finite increments Av;, AZ; and 
At 


’ 


Ag*® = d*At, 


where the superscripts k = 0,1,2,... indicate the progress of time, 
t = kAt. The position and velocity of a point 7 after time increment 
At are expressed as follows: 








ok+l _ =k ok 
U; = U;° + Ad; 

okt+l _ ok ~k 
x; 7 4 TT Av, 


The iterative computation of the velocities v,* and positions 7} is car- 


ried out for consecutive values of index k until all increments Av; and 
Az; fall below a threshold value. This indicates that the equilibrium 
state has been approximated to the desired accuracy. The next deriva- 
tion step is then performed. A system of equations corresponding to 
the resulting map is created, and the search for an equilibrium state 
resumes. In this way, the developmental process is simulated as periods 
of continuous cell expansion, delimited by instantaneous cell divisions. 
Continuity of cell shapes during divisions is preserved by the rule that 
sets the initial positions of vertices. 

For example, Figure 7.7 illustrates the expansion of a structure gen- 
erated by the L-system specified in Figure 7.3 in a derivation of length 
4. Figure 7.7a shows the structure immediately after the insertion of 
division walls, Figures 7.7b and c present intermediate wall positions, 
and Figure 7.7d describes the final structure at equilibrium. 


7.3. Microsorium linguaeforme 


In this section, the described simulation method is applied to visualize 
the development of the fern gametophyte Microsorium linguaeforme. 
The biological data is based on observations conducted by de Boer [22]. 
Fern gametophytes represent the sexually reproducing life stage of fern 


Apical activity 
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Figure 7.7: Expansion of a cellular structure: (a) after division, (b) and ( 
intermediate positions, (d) equilibrium 


plants. They show no differentiation into stem, leaf and root, forming 
a plant body called a thallus. The development of a thallus can be 
described conveniently in terms of the activity of the apical cell giving 
rise to segments, and the development of these segments. The modeling 
process captures repetitive patterns of cell divisions, so that large cel- 
lular structures can be described using a small number of productions. 

The apical cell is the originator of the gametophyte structure. It 
divides repetitively, giving rise each time to a new apical cell and a 
primary (initial) segment cell. The segment cells subsequently develop 
into multicellular segments. The division wall of an apical cell is at- 
tached to the thallus border on one side and to a previously created 
division wall on the other side. Thus, the division walls are oriented al- 
ternately to the left and to the right, yielding two columns of segments 
separated by a zig-zag dividing line (Figure 7.8). The recursive nature 
of the apical activity can be expressed by the following formula (called 
a cell division system [22]): 


A; — S| Ar Ar — Az |Sr 


This notation means that the cell on the left side of the arrow produces 
two daughter cells separated by a wall. 
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apical front C 


“ ” ” ” 
(0) (1) (2) (3) 


Figure 7.8: Apical production of segments. The labels Az and Az denote 
apical cells producing right segment Sp and left segment S;,, respectively. 
Dashed lines indicate the newly created division wall. The superscripts rep- 
resent segment age. The internal structure of segments is not shown. 
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Figure 7.9: Developmental sequence of a Microsorium segment 


In describing the structure of a segment, we distinguish between pericli- Division 
nal and anticlinal walls. Intuitively, periclinal walls are approximately pattern of 
parallel to the apical front of the thallus, and anticlinal walls are per- segments 


pendicular to this front. A more formal definition follows. 


e Ina primary segment, the apical front wall and one or more walls 
opposing it are periclinal walls. The remaining walls are anticlinal 
walls. 


e A division wall attached to two periclinal walls is an anticlinal 
wall, and vice-versa. 


In Microsorium, a wall is never attached to a periclinal wall on one 
side and an anticlinal wall on the other side, so the above definition 
comprises all possible cases. 

Microscopic observations of growing Microsorium gametophytes re- 
veal that most segments follow the same developmental sequence, shown 
diagrammatically in Figure 7.9. The primary segment cell Q, is first 
divided by a periclinal wall into two cells, Q2 and Q3. Subsequently, the 
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Figure 7.10: Developmental sequence of a Microsorium gametophyte 
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basal cell Q3 is divided by another periclinal wall into two “terminal” 
cells T that do not undergo further divisions. At the same time, the cell 
Qz lying on the thallus border is divided by an anticlinal wall into two 
cells of type Q;. Each of these cells divides in the same way as the pri- 
mary cell. Consequently, the recursive nature of segment development 
can be captured by the following cell division system: 


a2 Q2 > Q1| Q1 QF 


In the above rules, a horizontal bar denotes a periclinal wall between 
cells, and a vertical bar denotes an anticlinal wall. 

The development of the Microsorium thallus is a result of concur- 
rent divisions of the apical and segment cells. A single division of the 
apical cell corresponds to a single step in segment development. A de- 
velopmental sequence that combines the activity of the apex and the 
segments is shown in Figure 7.10. This figure also reveals offsets be- 
tween neighboring walls. It can be assumed that periclinal division 
walls form S-offsets in the segments on the right side of the apex and 
Z-offsets in the segments on the left side. 

In order to capture the development of Microsorium using the for- 
malism of map L-systems, it is necessary to identify all combinations 
of cells that may lie on both sides of a wall. Careful examination 
of these combinations yields the wall labeling scheme shown in Fig- 
ure 7.10. Two walls have the same label if and only if they divide in 
the same way.” The uppercase letters apply to right segment walls, and 
the corresponding lowercase letters denote symmetric walls in the left 
segments. A comparison of pairs of subsequent structures yields the 
L-system in the figure. 

Apical cell division results from the application of productions r, 
and Jl, (creation of a right segment) or J; and rz (creation of a left 
segment). Subsequent cell divisions in right and left segments proceed 
symmetrically. The development of a right segment is examined in 
detail. 

The insertion of wall segment B creates the first cell Q; of segment 
Sp (cf. Figure 7.8). Concurrently, wall D on the opposite side of the 
segment is transformed into F' (step 1). This transformation introduces 
a one-step delay before the application of production r4 which, together 
with rg, splits cell Q; into Q2 and Qs by the first periclinal wall H (step 
2). As the derivation progresses, production 4 inserts subsequent per- 
iclinal walls H between pairs of anticlinal walls F' (step 4). Production 





“Tt is conceivable to formulate an algorithm that would assign labels consistent 
with the above rule automatically. However, the labeling scheme given in Figure 7.10 
was obtained “by hand.” 
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Figure 7.11: Developmental sequence of a basal Microsorium segment 


r3 introduces a delay needed to create walls F’, which are inserted be- 
tween periclinal walls H and apical walls J using productions r; and rz. 
Production rg plays a role analogous to r3 — it introduces a one-step 
delay into the cycle which creates markers F at the apical front of the 
segment. Thus, periclinal walls H and anticlinal walls F’ are produced 
alternately, in subsequent derivation steps. The last two productions, 
rg and rg, create terminal walls x that do not undergo further changes. 
The first such wall is inserted between walls labeled G and E during 
derivation step 3. The subsequent walls x are inserted every second 
step between pairs of walls G; only production rg is applied in these 
cases. 

The L-system in Figure 7.10 was formulated under the assumption 
that all segments develop in the same way. However, in a real organism 
the two oldest segments, situated at the thallus base, form a modified 
pattern with less extensive cell divisions. The developmental sequence 
of a right basal segment is shown in Figure 7.11. The corresponding 
cell division system is given below. 


a2 Q2 > Qi | T QF 


The development of a Microsorium gametophyte including basal 
segments is captured by Figure 7.12. Only productions describing the 
development of the right side of the thallus are given. Their prede- 
cessors are denoted by uppercase letters. The corresponding lowercase 
productions, which complete the L-system definition, can be obtained 
by switching the “case” of letters and the orientation of markers. Wall 
directions remain unchanged. For example, the right-side production 
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Figure 7.12: The initial map and productions for the right side of a Microso- 
rium linguaeforme with basal segments 


corresponds to the left-side production 


be? Dis a[+Ble. 


Additional examples can be found by comparing the left and right 
columns of the L-system in Figure 7.10. 

Assuming the starting map specified by Figure 7.12, a simulated 
developmental sequence interpreted using the dynamic method to de- 
termine cell shape is given in Figure 7.13. Different colors are used 
to indicate the apical cell, the alternating “regular” segments, and the 
basal segments. A comparison of the last developmental stage with a 
photograph of a real gametophyte (Figure 7.14) shows good correspon- 
dence between the model and reality with respect to structure topology, 
the relative sizes and shapes of cells, and the overall shape of the thal- 
lus. This result is particularly interesting from a biological perspective, 
since it indicates that genetically controlled cell division patterns play 
an important role in determining the shape of a structure. 
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Figure 7.13: Simulated development of Microsorium linguaeforme 
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Figure 7.14: Microphotograph of Microsorium linguaeforme at magnification 
70x 
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Figure 7.15: Developmental sequence of a Dryopteris segment 


7.4 Dryopteris thelypteris 


Gametophytes of other fern species follow a similar developmental pat- 
tern, with the apex producing segments alternately to the left and 
right. However, the cell division patterns within segments vary between 
species yielding different thalli shapes. For example, Figure 7.15 shows 
segment development in Dryopteris thelypteris. The corresponding cell 
division system is given below. 


Qu D 
—> —— 


i= Qe O37 A Qs > Qi | Q1 


Qs Qs 


A developmental cycle of length 3, starting at cell Q1, produces 
two new cells Q; separated by an anticlinal wall and a sequence of 
four terminal cells T’ separated by periclinal walls. A developmental 
sequence that combines the activity of the apex and the segments is 
shown in Figure 7.17. As in the case of the map L-system in Figure 7.12, 
only productions for the right side of the thallus are shown. 

Apical cell division results from the application of productions r, 
and ly (creation of a right segment) or J; and r2 (creation of a left 
segment). The subsequent cell divisions proceed in a symmetric way 
in right and left segments. The development of a right segment is 
described below. 

The insertion of wall segment B creates the first cell Q; of segment 
Sp (cf. Figure 7.8). Concurrently, wall D on the opposite side of 
the segment is transformed by production rg into FG (step 1). This 
transformation introduces a one-step delay before the application of 
production r7 which, together with ro, splits cell Q; into Q2 and Q3 by 
the first periclinal wall x (step 2). Meanwhile, production r3 replaces 
wall C' by wall J (step 1), after which r4 replaces J by E (step 2). This 
introduces a two-step delay before cell Q3 is subdivided periclinally into 
two cells T by production r; (step 3). In the same step, productions r¢ 
and rg subdivide cell @2 into cells Q3 and Q4, separated by periclinal 
wall H. Wall O from step 1 is transformed into wall R by productions 
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Figure 7.16: Simulated development of Dryopteris thelypteris 


Ti5 (step 2) and rig (step 3). Walls R and H are replaced by ri7 and 
rja, resulting in the first anticlinal division of cell Q4 into two cells Q, 
by wall J (step 4). At the same time, productions r; and r¢ split cell 
Q3 periclinally into two cells 7. In the following derivation steps, each 
of the newly created cells Q, undergoes a sequence of changes similar 
to that described above. Production rg introduces a one-step delay 
before Q, is subdivided into Q2 and Qs3 using rg and rip (analogous 
to rg and rz). Productions rj; and rj play a role similar to r; and 
rg, while production 733 introduces a delay. Walls labeled x do not 
undergo further changes and cells T do not subdivide. A simulated 
developmental sequence generated by the L-system in Figure 7.17 using 
the dynamic method to determine cell shape is given in Figure 7.16. 

A comparison of Microsorium and Dryopteris gametophytes (either 
real or modeled) indicates that different division patterns of segment 
cells have a large impact on the overall thalli shapes. In Microsorium, 
the number of marginal cells, situated at the apical front of a segment, is 
doubled every second derivation step. The segments are approximately 
triangular, with a wide apical front, which results in the circular thallus 
shape. The apical front of Dryopteris segments is comparatively less 
developed. The number of marginal cells is doubled only every third 
step, and the segments grow faster in length. The resulting thallus 
shape is concave. 


Microsorium 
us. Dryopteris 


164 Chapter 7. Modeling of cellular layers 





(4) 


Figure 7.17: Developmental sequence of a Dryopteris gametophyte 
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Figure 7.17 (continued): Developmental 
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In order to quantify the relationship between segment cell division 
pattern and thallus shape, de Boer [22] proposed an empirical measure 
called periclinal ratio. It is based on the following considerations: 


e anticlinal growth takes place mainly along the margin and is ex- 
ponential since all marginal cells divide; and 


e periclinal growth of a segment is linear, as cells displaced away 
from the margin eventually stop dividing. 


Since the division pattern is recursive, the average ratio of the numbers 
of marginal cells in neighboring segments converges to a constant A. 
Similarly, the difference between the numbers of cells along the pericli- 
nal boundary of two neighboring segments converges to a constant P. 
By computer simulation [22], it was found that if the periclinal ratio 
P/A is smaller than 1.25, the apical front is convex as in the case of 
Microsorium linguaeforme. In contrast, if P/A is larger than 2.0, the 
simulated structure develops a concave apical front, corresponding to a 
heart-shaped thallus as in Dryopteris thelypteris. These results relate 
local growth to global shape. Periclinal ratios from 1.25 to 2.0 were not 
studied in detail. 


7.5 Modeling spherical cell layers 


Although the scope of this book is limited to plants, it is interesting 
to note that the formalism of L-systems can be applied also to simu- 
late some developmental processes in animals. For example, during the 
cleavage stage of development, an animal embryo consists of a single 
layer of cells that covers the surface of a spherical cavity. This structure 
is known as the blastula [6]. The cells divide synchronously in a regular 
pattern up to and including the 64-cell stage (6th cleavage). This devel- 
opment can be captured using an mBPMOL-system operating on the 
surface of a sphere rather than on a plane. To this end, cell walls are 
represented as great circle arcs connecting vertices that are constrained 
to the sphere surface. 

The extension of the dynamic interpretation method from the plane 
to the surface of a sphere requires few changes. Osmotic pressure and 
wall tension are calculated as before. Since the resulting force may 
displace a vertex away from the surface of the sphere, the actual vertex 
position is found by projecting the displaced point back to the sphere. 
During the cleavage stage, cells of embryos do not expand, thus the 
overall size of the sphere is constant. 

For example, deBoer [22] proposed the map L-system in Figure 7.18 
to model the development of a snail embryo, Patella vulgata, according 
to data presented by Biggelaar [150]. The starting map and develop- 
mental sequence are shown in Figure 7.18, while Figure 7.19 presents 
an alternative rendering. Each cell is approximated by a sphere cen- 
tered at the point obtained by raising the center of gravity of the cell 


7.5. Modeling spherical cell layers 167 





pr: A — bl[-alx|+alb 
D2 a — B[+A]x|—A]B 
P3 B-— a 

pe t> =k 

pp : © — D'I[+alE 

Pet D = C|=Alx 

P7 E = 

ps : F — E [—alG[+a] E 
po: GG J 

Dig? A = 1 

pu: I — B[—A]x[+4]B 
pi: J — bl[+a]x[—a]b 
pis: Z — C [—FIH[+F] c 


Figure 7.18: Developmental sequence of Patella vulgata (equatorial view) 
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vertices to the surface of the underlying spherical cavity. The radius is 
the maximum distance from the center of gravity to the cell vertices. 
A comparison of the Patella model at the 16-cell stage (bottom left 
of Figure 7.19) with an electron microscope image (Figure 7.20) shows 
good correspondence between the model and reality. 


7.6 Modeling 3D cellular structures 


The previous sections presented a method for modeling cellular lay- 
ers extending in a plane or on the surface of a sphere. However, real 
cellular structures are three-dimensional objects. In order to capture 
the three-dimensional aspect of cellular layers and model more complex 
structures, Lindenmayer [85] proposed an extension of map L-systems 
called cellwork L-systems. The notion of a cellwork is characterized as 
follows. 


e A cellwork is a finite set of cells. Each cell is surrounded by one 
or more walls (faces). 


e Each wall is surrounded by a boundary consisting of a finite, 
circular sequence of edges that meet at vertices. 


e Walls cannot intersect without forming an edge, although there 
can be walls without edges (in the case of cells shaped as spheres 
or tori). 


e Every wall is part of the boundary of a cell, and the set of walls 
is connected. 


e Each edge has one or two vertices associated with it. The edges 
cannot cross without forming a vertex, and there are no vertices 
without an associated edge. 


e Every edge is a part of the boundary of a wall, and the set of 
edges is connected. 


Note that the terms celland wall have different meanings in the context 
of cellworks than in the context of maps. 

The development of three-dimensional structures is captured using 
an extension of mBPMOL-systems called marker Binary Propagating 
Cellwork OL-systems |42]. An mBPCOL-system G is defined by a finite 
alphabet of edge labels 1, a finite alphabet of wall labels T, a starting 
cellwork w, and a finite set of edge productions P. The initial cellwork w 
is specified by a list of walls and their bounding edges. As in the case of 
mBPMOL-systems, edges may be directed or neutral. Each production 
is of the form A: 3 — a, where the edge A € © is the predecessor, the 
string 0 € {['*, x} is a list of applicable walls (* denotes all walls), and 
the string a is the successor, which is composed of edge labels from ©, 
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Figure 7.20: An electron microscope image of Patella vulgata 
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Figure 7.21: The phases of a derivation step 


wall labels from [ and the symbols | and |. The symbols outside square 
brackets describe the subdivision pattern for the predecessor. Pairs of 
matching brackets [ and | delimit markers that specify possible attach- 
ment sites for new edges and walls. As in the two-dimensional case, 
arrows indicate the relative directions of successor edges and markers 
with respect to the predecessor edge. The list @ contains all walls into 
which a marker should be inserted. In addition to the labels for edges 
and markers, a production successor specifies the labels of walls that 
may be created as a result of a derivation step. 

The syntax of a production is best explained using an example. The 
production 

A:14— DC2|Es]B3F 

applies to an edge A if it belongs to one or more walls labeled 1 or 
A (Figure 7.21a). The predecessor edge is subdivided into four edges 
D, C, B and F. During a derivation step, marker FE is introduced 
into all walls of type 1 or 4 that share edge A (Figure 7.21b), and 
can be connected with a matching marker inserted into the same wall 
by another production. As a result, the wall will split into two. The 
daughter wall having C’ as one of its edges will be labeled 2, and the 
wall having B as an edge will be labeled 3 (Figure 7.21c). Markers E 
can be connected only if both productions assign labels to the daughter 
walls in a consistent way. Otherwise, the markers are considered non- 
matching and are discarded. If several walls bounding a cell split in 
such a way that the sequence of new edges forms a closed contour, a 


7.6. Modeling 3D cellular structures 171 





Figure 7.22: Example of consistent edge productions 


new wall bound by these edges may be created. In order for this to 
occur, all markers involved must specify the same label for the new 
wall, 5 in this example (Figure 7.21d). 

The limitation of the scope of a production to specific walls may 
create a consistency problem while rewriting edges. For instance, as- 
sume that walls 1 and 2 share edge A and the following productions 
are in P: 


Pr: A:1 => CE 
pz: A:2 — AB 


Productions p; and pz are inconsistent since they specify two different 
partitions of the same edge. It is assumed that the mBPCOL-systems 
under consideration are free of such inconsistencies. This does not 
preclude the possibility of applying several productions simultaneously 
to the same edge. For example, a production pair, 


Diz A:l 3 Co[F'3] Es 
po: A:2 — Cs|DelEv, 


consistently divides edge A into segments C' and E, although the mark- 
ers inserted into walls 1 and 2 are different (Figure 7.22). 

According to the above discussion, a derivation step in an mBPCOL- 
system consists of three phases. 


e Each edge in the cellwork is replaced by successor edges and mark- 
ers using one or more productions in P. 


e Each wall is scanned for matching markers. If a match inducing 
a consistent labeling of daughter cells is found, the wall is subdi- 
vided. The selection of matching markers is done by the system. 
Unused markers are discarded. 


e Each cell is scanned for a circular sequence of new division edges. 
If a cycle assigning the same label to the division wall is found, it is 
used to bound the wall that will divide the cell into two daughter 
cells. If different possibilities exist, the edges are selected by the 
system. 
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(2) (3) 


Pi: A:1—-> By [A2| Bi 
Pp2 : A:2 > Bo [C2] Bo 
py: Bix* — A 


Figure 7.23: Example of a cellwork L-system 


A wall may be subdivided more than once as long as new division edges 
do not intersect and a consistent labeling of daughter walls is possible. 
In contrast, a cell may be divided only once in any derivation step. 

For example, Figure 7.23 presents a three-dimensional extension of 
the map L-system from Figure 7.2. In the first derivation step, produc- 
tion p; divides walls labeled 1, and production p2 divides walls labeled 
2. The inserted edges form a cycle that divides the cell with a new wall 
labeled 2. In the subsequent steps this process is repeated, generating 
a pattern of alternating division walls. Production p3 introduces the 
necessary delay. 

The dynamic method for interpreting map L-systems is extended to 
cellwork L-systems using the following assumptions: 


e the structure is represented as a three-dimensional network of 
masses corresponding to cell vertices, connected by springs which 
correspond to cell edges, 


e the springs are always straight and obey Hooke’s law, 


e for the purpose of force calculations, walls can be approximated 
by flat polygons, 


e the cells exert pressure on their bounding walls; the pressure on 
a wall is directly proportional to the wall area and inversely pro- 
portional to the cell volume, 
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e the pressure on a wall is divided evenly between the wall vertices, 
e the motion of masses is damped, and 
e other forces are not considered. 


The total force F; acting on a vertex X is given by the formula 


Fp= Sl F+ S- B+ Fi, 


ecE wEew 


where F, are forces contributed by the set of edges FE incident to X, 
Fy, are forces contributed by the set of walls W incident to X, and 
Fq = —bv is a damping force. The forces F, act along the cell edges and 
represent wall tension. The forces F,, are due to the pressure exerted 
by the cells on their bounding walls. The total force of pressure P 
exerted by a cell on a wall w has direction normal to w and is equal 
to p- A, where p is the internal cell pressure and A is the wall area. 
Calculation of the polygon area proceeds as in the two-dimensional 
case. The pressure p is assumed to be inversely proportional to the 
cell volume, p ~ V~!, which corresponds to the equation describing 
osmotic pressure (Section 7.2). The volume V of a cell is calculated by 
tesselating the cell into tetrahedra. The resulting differential equations 
are formed and solved as in the two-dimensional case. 

A division pattern that frequently occurs in epidermal cell structures 
is described by the L-system in Figure 7.24, based on a cyclic cellwork 
L-system (a slightly different formalism) proposed by Lindenmayer [85]. 
Productions pj, p2, pe and p7 are responsible for cell divisions, while the 
remaining productions introduce delays such that the division pattern 
is staggered. 

On the surface, the cellular structures analyzed in this chapter may 
appear quite unrelated to the models discussed previously. However, 
a closer inspection reveals many analogies. For example, consecutively 
created segments of a fern gametophyte exhibit a phase effect corre- 
sponding to that observed in inflorescences. Furthermore, parts of an 
older gametophyte situated near the apex have the same topology as 
the entire thallus at an earlier developmental stage, which associates 
the recursive structures generated by map or cellwork L-systems with 
self-similar patterns created using string L-systems. As observed by 
Oppenheimer [105], self-similarity appears to be one of the general prin- 
ciples organizing the world of botany. The next chapter discusses this 
topic in more detail. 


Epidermal cells 
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Pi: A:123 —- C3[E1|Bo[D4|C3 
p2: A:4 = CBa|Fi|Ca 

p3: Bix — A 

pa: Ci — B 

py: Ei — D 

po: F:123 — HGH 

P7: F:4 — Ha|F1|Ga[F1]Ha 
pg: Gi — F 

po: H: * — G 


Figure 7.24: Developmental sequence of epidermal cells: (a) The starting 
cellwork; (b), (d) and (f) cellworks immediately after cell divisions; (c), (e) 
and (g) the corresponding cellworks at equilibrium 


Chapter 8 


Fractal properties of plants 


What is a fractal? In his 1982 book, Mandelbrot defines it as a set with 
Hausdorff-Besicovitch dimension Dy strictly exceeding the topological 
dimension Dr [95, page 15]. In this sense, none of the figures presented 
in this book are fractals, since they all consist of a finite number of 
primitives (lines or polygons), and Dy = Dr. However, the situation 
changes dramatically if the term “fractal” is used in a broader sense [95, 
page 39]: 


Strictly speaking, the triangle, the Star of David, and the 
finite Koch teragons are of dimension 1. However, both 
intuitively and from the pragmatic point of view of the sim- 
plicity and naturalness of the corrective terms required, it is 
reasonable to consider an advanced Koch teragon as being 
closer to a curve of dimension log 4/log 3 than to a curve 
of dimension 1. 


Thus, a finite curve can be considered an approximate rendering 
of an infinite fractal as long as the interesting properties of both are 
closely related. In the case of plant models, this distinctive feature is 
self-similarity. 

The use of approximate figures to illustrate abstract concepts has a 
long tradition in geometry. After all, even the primitives of Euclidean 
geometry — a point and a line — cannot be drawn exactly. An in- 
teresting question, however, concerns the relationship between fractals 
and real biological structures. The latter consist of a finite number of 
cells, thus are not fractals in the strict sense of the word. To consider 
real plants as approximations of “perfect” fractal structures would be 
acceptable only if we assumed Plato’s view of the supremacy of ideas 
over their mundane realization. A viable approach is the opposite one, 
to consider fractals as abstract descriptions of the real structures. At 
first sight, this concept may seem strange. What can be gained by 
reducing an irregular contour of a compound leaf to an even more ir- 
regular fractal? Would it not be simpler to characterize the leaf using 
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a smooth curve? The key to the answer lies in the meaning of the term 
“simple.” A smooth curve may seem intuitively simpler than a fractal, 
but as a matter of fact, the reverse is often true [95, page 41]. Accord- 
ing to Kolmogorov [|80], the complexity of an object can be measured 
by the length of the shortest algorithm that generates it. In this sense, 
many fractals are particularly simple objects. 

The above discussion of the relationship between fractals and plants 
did not emerge in a vacuum. Mandelbrot [95] gives examples of the re- 
cursive branching structures of trees and flowers, analyzes their 
Hausdorff-Besicovitch dimension and writes inconclusively “trees may 
be called fractals in part.” Smith [136] recognizes similarities between 
algorithms yielding Koch curves and branching plant-like structures, 
but does not qualify plant models as fractals. These structures are pro- 
duced in a finite number of steps and consist of a finite number of line 
segments, while the “notion of fractal is defined only in the limit.” Op- 
penheimer [105] uses the term “fractal” more freely, exchanging it with 
self-similarity, and comments: “The geometric notion of self-similarity 
became a paradigm for structure in the natural world. Nowhere is this 
principle more evident than in the world of botany.” The approach pre- 
sented in this chapter, which considers fractals as simplified abstract 
representations of real plant structures, seems to reconcile these previ- 
ous opinions. 

But why are we concerned with this problem at all? Does the no- 
tion of fractals provide any real assistance in the analysis and modeling 
of real botanical structures? On the conceptual level, the distinctive 
feature of the fractal approach to plant analysis is the emphasis on 
self-similarity. It offers a key to the understanding of complex-looking, 
compound structures, and suggests the recursive developmental mech- 
anisms through which these structures could have been created. The 
reference to similarities in living structures plays a role analogous to the 
reference to symmetry in physics, where a strong link between conser- 
vation laws and the invariance under various symmetry operations can 
be observed. Weyl [159, page 145] advocates the search for symmetry 
as a cognitive tool: 


Whenever you have to deal with a structure-endowed entity 
yi, try to determine its group of automorphisms, the group 
of those element-wise transformations which leave all struc- 
tural relations undisturbed. You can expect to gain a deep 
insight into the constitution of © in this way. 


The relationship between symmetry and self-similarity is discussed 
in Section 8.1. Technically, the recognition of self-similar features of 
plant structures makes it possible to render them using algorithms de- 
veloped for fractals as discussed in Section 8.2. 
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T,(A) 


Figure 8.1: The Sierpijski gasket is closed with respect to transformations 
T,, Tz and T3 (a), but it is not closed with respect to the set including the 
inverse transformations (b). 


8.1 Symmetry and self-similarity 


The notion of symmetry is generally defined as the invariance of a con- 
figuration of elements under a group of automorphic transformations. 
Commonly considered transformations are congruences, which can be 
obtained by composing rotations, reflections and translations. Could 
we extend this list of transformations to similarities, and consider self- 
similarity as a special case of symmetry involving scaling operations? 

On the surface, this seems possible. For example, Wey] {159, page 68] 
suggests: “In dealing with potentially infinite patterns like band orna- 
ments or with infinite groups, the operation under which a pattern is 
invariant is not of necessity a congruence but could be a similarity.” 
The spiral shapes of the shells Turritella duplicata and Nautilus are 
given as examples. However, all similarities involved have the same 
fixed point. The situation changes dramatically when similarities with 
different fixed points are considered. For example, the Sierpinski gasket 
is mapped onto itself by a set of three contractions T;, T> and T3 (Fig- 
ure 8.la). Each contraction takes the entire figure into one of its three 
main components. Thus, if A is an arbitrary point of the gasket, and 
T =T;,,T;,...T;,, is an arbitrary composition of transformations T), T> 
and 73, the image T(A) will belong to the set A. On the other hand, 
if the inverses of transformations 7), T> and 73 can also be included 
in the composition, one obtains points that do not belong to the set A 
nor its infinite extension (Figure 8.1b). This indicates that the set of 
transformations that maps A into itself forms a semigroup generated 
by T,, T> and T3, but does not form a group. Thus, self-similarity is a 
weaker property than symmetry, yet it still provides a valuable insight 
into the relationships between the elements of a structure. 
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Figure 8.2: The fern leaf from Barnsley’s model [7] 


8.2 Plant models and iterated function sys- 
tems 


Barnsley [7, pages 101-104] presents a model of a fern leaf (Figure 8.2), 
generated using an iterated function system, or IFS. This raises a ques- 
tion regarding the relationship between developmental plant models 
expressed using L-systems and plant-like structures captured by IF Ses. 
This section briefly describes IFSes and introduces a method for con- 
structing those which approximate structures generated by a certain 
type of parametric L-system. The restrictions of this method are ana- 
lyzed, shedding light on the role of IFSes in the modeling of biological 
structures. 

By definition [74], a planar iterated function system is a finite set 
of contractive affine mappings T = {7\,7>,...,7;,} which map the 
plane R x R into itself. The set defined by T is the smallest nonempty 
set A, closed in the topological sense, such that the image y of any 
point « € A under any of the mappings J; € JT also belongs to A. 
It can be shown that such a set always exists and is unique [74] (see 
also [118] for an elementary presentation of the proof). Thus, starting 
from an arbitrary point x € A, one can approximate A as a set of 
images of x under compositions of the transformations from 7. On 
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Figure 8.3: A comparison of three attracting methods for the rendering of 
a set defined by an IFS: (a) deterministic method using a balanced tree of 
depth n = 9 with the total number of points Ny = 349, 525, (b) deterministic 
method using a non-balanced tree with N2 = 198,541 points, (c) stochastic 
method with N3 = No points 
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the other hand, if the starting point x does not belong to A, the con- 
secutive images of x gradually approach A, since all mappings 7; are 
contractions. For this reason, the set A is called the attractor of the 
IFS 7. The methods for rendering it are based on finding the images 
Die Tie oie (TG(b)) cee) ) Bd pans digg Ly, ond “are termed atiracting 
methods. According to the deterministic approach [123], a tree of trans- 
formations is constructed, with each node representing a point in A. 
Various strategies, such as breadth-first or depth-first, can be devised to 
traverse this tree and produce different intermediate results [60]. If the 
transformations in J do not have the same scaling factors (Lipschitz 
constants), the use of a balanced tree yields a non-uniform distribu- 
tion of points in A. This effect can be eliminated by constructing a 
non-balanced tree, using a proper criterion for stopping the extension 
of a branch [60]. An alternative approach for approximating the set A 
is termed the chaos game [7| (see also [107, Chapter 5]). In this case, 
only one sequence of transformations is constructed, corresponding to 
a single path in the potentially infinite tree of transformations. The 
transformation applied in each derivation step is selected at random. 
In order to achieve a uniform distribution of points in the attractor, 
the probability of choosing transformation T; € T is set according to 
its Lipschitz constant. Figure 8.3 illustrates the difference between the 
stochastic and deterministic methods of rendering the attractor. The 
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underlying IFS consists of four transformations, given below using ho- 
mogeneous coordinates [40]: 


0.00 0.00 0.0 
0.00 0.16 0.0 
0.00 0.00 1.0 


ieee 
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0.20 0.23 
—0.26 0.22 
0.00 1.60 


—0.15 0.26 
0.28 0.24 
0.00 0.44 


0.85 —0.04 
0.04 0.85 
0.00  =1.60 


Th = 





T3 = 





Ty = 
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Other methods for the rendering of the set A, defined by an inter- 
ated function system 7, include the repelling or escape-time method 
and the distance method [60, 118]. Both methods assign values to 
points outside of A. The first method determines how fast a point is 
repelled from A to infinity by the set of inverse transformations T,', 
where 7; € J. An example of the application of this method, with 
escape time values represented as a height field, is shown in Figure 8.4. 
The second method computes the Euclidean distance of a point from 
the attractor A. 

The problem of constructing an IFS that will approximate a branch- 
ing structure modeled using an L-system can now be considered. This 
discussion focuses specifically on structures that develop in a biologi- 
cally justifiable way, by subapical branching (Section 3.2). The com- 
pound leaf shown in Figure 5.11a on page 129 will be used as a working 
example. In this case, the apical delay D is equal to zero, and the 
L-system can be represented in the simplified form: 


wi A 
py: A > * — F(1)[+A][-A]F(1)A (8.1) 
po: F(a) : * > F(axR) 


Figure 8.4: Fern dune —> 
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Figure 8.5: Initial sequences of structures generated by the L-systems spec- 
ified in equations (8.1) and (8.2) 


This L-system operates by creating segments of constant size, then 
increasing their length by constant factor R in each derivation step 
(Figure 8.5a). As discussed in Section 1.10.3, a structure with the 
same proportions can be obtained by successively appending segments 
of decreasing length (Figure 8.5b): 


w: A(1) (8.2) 
pi: A(s) : * — F(s)[+A(s/R)][—A(s/R)]F(s)A(s/R) 


Let A,,(s) denote the structure generated by module A(s) inn > 1 
derivation steps. According to production p;, the following equality 
holds: 


An(8) = F(s)[+An-1(8/)][—An-1(8/P)]F(s)An-1(s/R) (8.3) 
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Figure 8.6: Illustration of equation (8.4), with uw, = F(1) 
[—(45) F'(0.5)| F'(0.5) 
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It is important to clearly distinguish between a parametric word p (a 
string of modules) and its turtle interpretation 7(j) (a set of points in 
the plane). The symbol M(,.) will be used to denote the transformation 
induced by pw. This transformation moves the turtle from its initial 
position and orientation to those resulting from the interpretation of 
word yu. According to the definition of turtle interpretation (Chapter 1), 
if a word yp is decomposed into subwords 4, and pz such that 2 does 
not contain unbalanced right brackets, then 


T(t) = FT (usta) = TF (ba) UT (2) M (qa). (8.4) 


See Figure 8.6 for an illustration. By applying equation (8.4) to (8.3), 
we obtain 


T(An(8)) = 


S/RM 
which is true for any n > 1. Now let A(s) be the limit of the sequence 
of sets 7(A,(s)) from equation 8.6, 


A(s) = lim J(A,(s)). 


n—co 


At the limit we obtain: 


A(s) = J(F(2s)) U A(s/R\(Ti UTLU TS), (8.6) 
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where Tj = T(F(s)+), T5 = T(F(s)-), and T; = T(F(2s)). Let 
S(s/R) be the operation of scaling by s/R, then 


A(s/R) = A(s)S(s/R). 


By noting T; = T/S(s/R) for i = 1,2,3, equation (8.6) can be trans- 
formed to 


A(s) = J(F(2s)) U A(s)(Ty U Tp U Th). (8.7) 


The solution of this equation with respect to A(s) is 


where (7; U 72 U T3)* stands for the iteration of the union of transfor- 
mations T\, 7 and 73. Equation (8.8) suggests the following method 
for constructing the set (A(s)): 


e create segment 7(F'(2s)) 


e create images of 7(f(2s)) using transformations 7), 7>,73 and 
their compositions 


Equation (8.7) and the method of constructing the set A(s) based 
on equation (8.8) are closely related to the definition of iterated function 
systems stated at the beginning of this chapter. However, instead of 
starting from an arbitrary point x € A(s), the iteration begins with 
the set 7(F(2s)). Although this is simply a straight line segment, a 
question arises as to how its generation can be incorporated into an 
IFS. Two approaches can be distinguished. 

The first approach is related to the notions of hierarchical iterated 
function systems discussed by Reuter [123] and recurrent IF Ses intro- 
duced recently by Barnsley [8]. The line segment 7(F'(2s)) is gener- 
ated using an IFS, for example consisting of two scaling transformations 
Q, and Q» which map it onto its upper and lower half (Figure 8.7a). 
Subsequently, transformations 7), T> and 73 are applied to create other 
points of the set A(s). The order of transformation application is im- 
portant. Transformations Q; and Q2 are used solely for the purpose 
of initial segment creation. They must not be applied after 7\, 7) or 
T3, since in this case they would affect the branching structure under 
consideration. The admissible sequences of transformations can be de- 
fined using a directed control graph (Figure 8.7b), and correspond to 
the infinite set of paths starting at node a.! The term controlled it- 
erated function system (CIFS) denotes an IFS with restrictions on the 
transformation sequences imposed by a control graph. Thus, noting 
the angle increment associated with symbols + and — by 0, the fractal 





‘Formally, the sequence of admissible transformations is the regular language 
accepted by the finite (Rabin-Scott) automaton represented by the graph in Fig- 
ure 8.7b. 
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Figure 8.7: Construction of the set A(S): (a) definition of an IFS {Q1, Qa} 
that generates the initial line segment, (b) the control graph specifying the 
admissible sequences of transformation application 


a 


approximation of the leaf in Figure 5.11a is given by the CIFS with the Resulting CIFS 
control graph in Figure 8.7b and the transformations specified below: 
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The second approach to the generation of the line segment 7(F'(2s)) Noninvertible 
is consistent with the method applied by Barnsley to specify the fern transforma- 
leaf in Figure 8.2. The idea is to map the entire branching structure tions 


A(s) onto the line 7(F'(2s)). This can be achieved using a noninvertible 
transformation Q which collapses all branches into a vertical line. The 
scaling factor along the y axis is the ratio of the desired segment length 
2s, and the limit height of the entire structure A(s), 
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a b 











Figure 8.8: Two renderings of the compound leaf from Figure 5.11a, gener- 
ated using iterated function systems 


This last value is calculated as the limit of the geometric series with the 
first term equal to 2s and the ratio equal to 1/R. Thus, the compound 
leaf of Figure 5.11a is defined by an IFS consisting of transformation 


0 oO 0 
Q=|0 1-1/R 0 
0 oO 1 





and transformations 7\, T> and 73 specified as in the case of the con- 
trolled IFS. 

Two fractal-based renderings of the set A(s) are shown in Fig- 
ure 8.8. Figure 8.8a was obtained using the controlled IFS and a deter- 
ministic algorithm to traverse the tree of admissible transformations. 
Figure 8.8b was obtained using the “ordinary” IFS and the random 
selection of transformations. Figure 8.9 shows another fractal-based 
rendering of the same structure. The spheres have radii equal to the 
distance from the sphere center to the leaf, within a specified e. 


Figure 8.9: Carrot leaf — 
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An equation expressing the limit set as a union of 
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An equation expressing the limit object as the im- 
age of an initial object under an iteration of a 
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A (controlled) IFS 











Figure 8.10: Steps in the construction of an IFS given an L-system capturing 
a developmental model 


8.2. Plant models and iterated function systems 189 


It is instructive to retrace the logical construction that started with 
an L-system, and ended with an iterated function system which can 
generate fractal approximations of the same object (Figure 8.10). An 
analysis of the operations performed in the subsequent steps of this 
construction reveals its limitations, and clarifies the relationship be- 
tween strictly self-similar structures and real plants. The critical step 
is the transformation of the L-system with elongating internodes to the 
L-system with decreasing apices. It can be performed as indicated in 
the example if the plant maintains constant branching angles as well as 
fixed proportions between the mother and daughter segments, indepen- 
dent of branch order. This, in turn, can be achieved if all segments in 
the modeled plant elongate exponentially over time. These are strong 
assumptions, and may be satisfied to different degrees in real plants. 
Strict self-similarity is an abstraction that captures the essential prop- 
erties of many plant structures and represents a useful point of reference 
when describing them in detail. 


Conclusions 





Epilogue 


This quiet place, reminiscient of Claude Monet’s 1899 painting Water- 
lies pool — Harmony in green, does not really exist. The scene was 
modeled using L-systems that captured the development of trees and 
water plants, and illuminated by simulated sunlight. It is difficult not 
to appreciate how far the theory of L-systems and the entire field of 
computer graphics have developed since their beginnings in the 1960’s, 
making such images possible. Yet the results contained in this book are 
not conclusive and constitute only an introduction to the research on 
plant modeling for biological and graphics purposes. The algorithmic 
beauty of plants is open to further exploration. 


+ Figure E.1: Water-lilies 


Appendix A 


Software environment for 
plant modeling 


This book is illustrated with images of plants which exist only as math- 
ematical models visualized by means of computer graphics. The soft- 
ware environment used to construct and experiment with these models 
includes dozens of programs and hundreds of data files. This creates 
the nontrivial problem of organizing all components for easy definition, 
saving, retrieval and modification of the models. In order to solve it, the 
idea of simulation was extended beyond the level of individual plants 
to an entire laboratory in botany [98]. Thus, a user can create and 
conduct experiments in a virtual laboratory by applying intuitive con- 
cepts and techniques from the “real” world. As an operating system 
defines the way a user perceives a computing environment, the virtual 
laboratory determines a user’s perception of the environment in which 
simulated experiments take place. In the future, a virtual laboratory 
may complement, extend, or even replace books as a means for gather- 
ing and presenting scientific information. Because of this potential, the 
laboratory in which the research reported in this book was produced is 
described here in more detail. 


A.1 A virtual laboratory in botany 


A virtual laboratory, like its “real” counterpart, is a playground for ex- 
perimentation. It comes with a set of objects pertinent to its scientific 
domain (in this case, plant models), tools which operate on these ob- 
jects, a reference book and a notebook. Once the concepts and tools are 
understood, the user can expand the laboratory by adding new objects, 
creating new experiments, and recording descriptions in the notebook. 
An experienced user can expand the laboratory further by creating and 
installing new tools. 


User’s 
perspective 


Laboratory = 
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Requirements 


194 Appendix A. Software environment for plant modeling 


Technically, a virtual laboratory is a microworld which can be ex- 
plored under the guidance of a hypertext system. The term “microworld” 
denotes an interactive environment for creating and conducting simu- 
lated experiments. The guidance could be provided in the form of a 
traditional book, but an electronic document is more suitable for inte- 
gration with a microworld. In a sense, both components of the virtual 
laboratory are described by Nelson in Dream Machines |104]. The pi- 
oneering role of this book in introducing the concept of hypertext is 
known, but under the heading The Mind’s Eye the notion of a mi- 
croworld is also anticipated: 


Suppose that you have a computer. 
What sorts of things would you do with it? 
Things that are imaginative 
and don’t require too much else. 
I am hinting at something. 
You could have it make pictures and show you stuff 
and change what it shows depending on what you do. 


A virtual laboratory can be divided into two components: the ap- 
plication programs, data files and textual descriptions that comprise 
the experiments; and the system support that provides the framework 
on which these domain-dependent experiments are built. The following 
list specifies the features of this framework. 


e Consistent organization of the lab. In the lab environment, 
experiments are run by applying tools (programs) to objects (data 
files). An object consists of files that are grouped together so 
that they can be retrieved easily. The format of the objects is 
sufficiently standardized to allow straightforward implementation 
of common operations such as object saving and deletion. 


e Inheritance of features. It is often the case that several ob- 
jects differ only in details. For example, two lilac inflorescences 
may differ only in the color of their petals. The mechanism of 
inheritance is employed to store such objects efficiently. 


e Version control. Interaction with an object during experimen- 
tation may result in a temporary or permanent modification. In 
the latter case, the user is able to decide whether the newly cre- 
ated object replaces the old one or should be stored as another 
version of the original object. 


e Interactive manipulation of objects. The laboratory pro- 
vides a set of general-purpose tools for manipulating object pa- 
rameters. For example, objects can be modified using control 
panels or by editing specific fields in a textual description of an 
experiment. 
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e Flexibility in conducting experiments. The user may apply 
tools to objects in a dynamic way while an experiment is being 
conducted. This can be contrasted to a static experiment de- 
signed when the object is initially incorporated into the system. 


e Guidance through the laboratory. A hypertext system im- 
poses a logical organization on the set of objects, provides a 
textual description of the experiments, and makes it possible to 
browse through the experiments in many ways. Specific experi- 
ments are invoked automatically when the corresponding text is 
selected, in order to facilitate demonstrations and assist a novice 
user. 


So far, objects have been referred to in an intuitive way, relying 
on the analogy between a real and virtual laboratory. For example, if 
our interest is in the development of the gametophyte Microsorium lin- 
guaeforme, in a real laboratory we would experiment with a specimen 
of the plant, while in a virtual laboratory we explore the corresponding 
mathematical model. However, the analogy to real objects does not ex- 
tend to the level of detailed object definition. Specific design decisions 
are needed for software development purposes. In the current design, a 
laboratory object is defined as a directory containing two types of files 
and a subdirectory. 


e The data files comprise our knowledge of a particular model. 


e A specification file defines the data files which make up the object 
and the tools which apply to them. 


e A directory of extensions lists objects which inherit some features 
of the current object. 


The object-oriented file structure which provides the basis for lab op- 
eration can be represented by a hierarchy of directories and files (Fig- 
ure A.1). 

The path of subdirectories leading to an object establishes the in- 
heritance structure for the lab. Inheritance is based on the idea of 
specifying new objects in reference to objects which already exist [81]. 
The “old” object is called a prototype and the new one is its extension. 
The extension contains only those files which are different from the 
corresponding files in the prototype. Files that remain the same are 
delegated to the prototype by establishing links. In other words, the 
object directory will contain those files that are unique to the object, 
and links to files that are inherited from its prototype (Figure A.2). 
This approach saves space, facilitates creation of objects similar to the 
prototype, and allows a single change in the prototype to propagate 
through all descendents. 
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Figure A.1: The hierarchical structure of objects 
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Figure A.2: A prototype and its extension. Shaded areas indicate links. 
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Figure A.3: An object icon with menus 














To conduct an experiment, all files that make up the selected object 
are copied to a temporary location called the lab table. Consequently, 
manipulation of object parameters does not disturb the stored version. 
When the experiment is finished, the user may save the results by 
overwriting the original object or by creating an extension. In the 
latter case, the files on the lab table are compared with those in the 
prototype object; those files that differ from the prototype are saved, 
and links to the remaining files are established automatically. 

The ability to manipulate the parameters in an experiment easily is 
an essential feature of the virtual laboratory. As a rule, all parameters 
involved in an experiment are supplied to the tools through the object’s 
data files. In order to modify a parameter, the user edits the appropri- 
ate file, which is subsequently re-read by the application. Though the 
editing of parameters can be accomplished using a text editor, in many 
cases parameter modification can be performed more conveniently using 
virtual control panels [114]. The current implementation of the labora- 
tory provides the user with a general-purpose control manager which 
creates panels according to user-supplied configuration files. 

The user is able to apply a tool to an object as a whole, without 
detailed knowledge of the programs involved or the component files. 
This is achieved through the object’s specification file which lists all 
files associated with an object and the tools that can be applied to 
them. This information is used to create a hierarchy of menus associ- 
ated with an icon representing the object (Figure A.3). The end nodes 
in the hierarchy invoke tools that operate on the object. For example, 
selection of the item image followed by the item generate from the 
menus in the figure would invoke the plant modeling program Pfg. 

A user may browse through the objects in the lab by following either 
the hierarchical structure of objects or hypertext links. The browser is 
used to navigate through the hierarchy, moving down through succes- 
sive extensions or up through previous levels. At any time, the user 
may request that an object be placed on the lab table. The hypertext 
document associated with the lab provides an alternative method of 
browsing and a means of relating objects independent of the hierarchy. 


Version control 


Object 
manipulation 


Tool 
application 


Browsing 


198 Appendix A. Software environment for plant modeling 


A.2 List of laboratory programs 


The essential programs incorporated into the virtual laboratory in botany 
are listed below. 


e Plant and fractal generator (Pfg) 

P. Prusinkiewicz and J. Hanan 

Given an L-system, a set of viewing parameters and optional files 
specifying predefined surfaces, Pfg generates the modeled struc- 
ture by carrying out the derivation, then interpreting the resulting 
string using turtle geometry. Both non-parametric and paramet- 
ric L-systems are supported. The model can be visualized directly 
on the screen of an IRIS workstation or output to a file. The first 
mode of operation is used to experiment with the model interac- 
tively and present developmental sequences. The output file can 
be either in Postscript format, particularly suitable for printing 
results such as fractal curves and inflorescence diagrams on a laser 
printer, or in the format required by the ray-tracer Rayshade for 
realistic rendering of the modeled structures. 


e Modeling program for phyllotactic patterns (Spiral) 

D. R. Fowler 

Spiral is an interactive program for modeling organs with spiral 
phyllotactic patterns. The user can choose between planar and 
cylindrical patterns, and modify parameters which define model 
geometry (Chapter 4). This technique is faster than “growing” 
organs using parametric L-systems. Once an organ has been de- 
signed, it can be expressed using an L-system and incorporated 
into a plant structure. 


e Interactive surface editor (Ise) 
J. Hanan 
Ise makes it possible to define and modify bicubic surfaces con- 
sisting of one or several arbitrarily connected patches. The output 
files produced by Ise are compatible with Pfg and Spiral. 


e Modeling program for cellular structures (Mapl) 

F. D. Fracchia 

Mapl accepts the specification of a two-dimensional cell layer cap- 
tured by a map L-system and generates the resulting developmen- 
tal sequence using the dynamic method of map interpretation. 
Options include map generation on the surface of a sphere, and 
the simulation of development in three dimensions according to 
a given cellwork L-system. As in the case of Pfg, the models can 
be visualized directly on the screen or output to a file in either 
Postscript or Rayshade format. 
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e Control panel manager (Panel) 

L. Mercer and A. Snider 

This program creates control panels containing sliders and but- 
tons, according to a configuration file provided by the user. Upon 
activation of a control by the mouse, Panel generates a message 
which indicates the corresponding control value. Application pro- 
grams process this information and modify the appropriate pa- 
rameters. For example, a panel can be used to control parameters 
used by Pfg, Spiral, or Mapl. 


e Ray tracer (Rayshade) 
C. Kolb, Yale University 


Rayshade reads a scene description from a text file, and renders it 
using ray tracing. Scenes can be composed of primitives such as 
planes, triangles, polygons, spheres, cylinders, cones and height 
fields, grouped together to form objects. These objects can be 
instantiated in other object definitions to create a hierarchical 
description of a scene. Transformations including translation, ro- 
tation and scaling, and a variety of procedural textures can be ap- 
plied to any object. Extended light sources, simulation of depth 
of field, and adaptive supersampling are supported. The program 
uses 3D grids to partition object space for fast intersection tests. 


e Previewer for the ray tracer (Preray) 
A. Snider 
Preray is a previewer for Rayshade used to provide a fast wire 
frame rendering of a scene before committing time to ray tracing. 
A control panel associated with Preray makes it possible to set 
viewing parameters interactively. 


e Modeling program based on Euclidean constructions (L.E.G.O.) 
N. Fuller 
L.E.G.O. makes it possible to model two- and three-dimensional 
objects using geometric constructions. In the scope of this book, 
L.E.G.O. was used to model man-made objects such as the Zinnia 
vase and the Water-lilies bridge. 


e Iterated function system generator (Ifsg) 

D. Hepting 

A fractal defined by an iterated function system is described by 
a finite set of contractive affine transformations with an optional 
finite state control mechanism. /fsg accepts input from a file spec- 
ifying the transformations and rendering information. The pro- 
gram is capable of rendering by either attracting, distance-based 
or escape-time methods. The output can be displayed directly 
on an IRIS workstation or written to a file for further processing. 
In the scope of this book, [fsg was used to obtain results which 
related plant models expressed using L-systems to fractals. 
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Figure A.4: A virtual laboratory screen 


Figure A.4 presents a sample screen of a Silicon Graphics IRIS 4D/60 
workstation running some of the above programs within the virtual 
laboratory framework. The icon in the top right corner represents the 
laboratory browser which was used to select a sunflower plant as the 
current object. The icon underneath and the associated menu were 
subsequently applied to select tools which operate on the object. The 
control panel in the bottom right corner of the screen is a part of the 
surface editor Ise. The manipulated petal is displayed as a wire frame 
in the window labeled Jse, and incorporated into a flower head by the 
modeling program Spiral which presents its output in the window sun- 
flower. The flower heads are in turn incorporated into a complete plant 
model generated by Pfg and rendered using Rayshade in the window 
plant.rle. The panel below that window makes it possible to choose 
organs included in the model and change parameters related to the 
angles of the branching structure. The metaphor of a virtual labora- 
tory provides a uniform interface to various operations on the selected 
plant, ranging from the modification of a petal to the rendering of the 
complete model. 


Appendix B 


About the figures 


The following descriptions of the color images include details about the 
pictures not described in the main text. Unless otherwise stated, figures 
were created at the University of Regina. 


Figure 1.19 [page 20] Three-dimensional Hilbert curve 
F. D. Fracchia, P. Prusinkiewicz, N. Fuller 
(1989) 


This image was rendered using ray-tracing without shadows. 


Figure 1.25 [page 26] Three-dimensional bush 
P. Prusinkiewicz (1986) 


Simple branching structure, rendered using the firmware of a Sil- 
icon Graphics IRIS workstation. Total generating and rendering 
time on IRIS 4D/20: 4 seconds. 


Figure 1.28 [page 29] Flower field 
P. Prusinkiewicz (1986) 


The field contains four rows of four plants. The scene was ren- 
dered with IRIS firmware, using depth-cueing to assign colors to 
petals. 


Figure 1.35 [page 45] Developmental stages of Anabaena catenula 
J. Hanan, P. Prusinkiewicz (1989) 


Figure 2.1 [page 52] Organic architecture 
Ned Greene, NYIT (1989) 
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An array of 300 x 300 x 300 voxel space automata was used to 
track a polygonal model of a house. Rendering was performed 
using a probabilistic radiosity method. See [54] for a full descrip- 
tion. 


Figure 2.3 [page 54] Acer graphics 


Jules Bloomenthal, NYIT (1984) 


A model of a maple tree. The basic branching structure was gen- 
erated recursively. Limbs were modeled as generalized cylinders, 
obtained by moving discs of varying radii along spline curves. 
Real bark texture was digitized and used as a bump map. Leaf 
texture was obtained by digitizing a photograph of a real leaf and 
emphasizing the veins using a paint program. See [11] for details. 


Figure 2.4 [page 54] Forest scene 


Bill Reeves, Pixar (1984) 


A scene from the film The Adventures of André and Wally B, 
modeled using particle systems. Shading and shadows were ap- 
proximated using probabilistic techniques. Visible surfaces were 
determined using depth-sorting. See [119] for a full description. 


Figure 2.5 [page 55] Oil palm tree canopy 


CIRAD Modelisation Laboratory (1990) 


A developmental model of oil palm trees, modeled using the 
method originated by de Reffye and described from the graph- 
ics perspective in [30]. 


Figure 2.10 [page 61] Medicine lake 


F. K. Musgrave, C. E. Kolb, P. Prusinkiewicz, 
B. B. Mandelbrot (1988) 


A scene combining a fractal terrain model, a tree generated us- 
ing L-systems, and a rainbow. The rainbow model was derived 
from a simulation of refraction with dispersion of light through 
an idealized raindrop. Procedural textures were applied to the 
mountains, the water surface and a vertical plane modeling the 
sky. See [101] for further details. 


Figure 2.11 [page 62] Surrealistic elevator 


A. Snider, P. Prusinkiewicz, N. Fuller (1989) 


The elevator was modeled using L.E.G.O. The island is a su- 
perquadratic surface. Procedural textures were applied to create 
stars in the sky, craters on the moon, colored layers in the rock, 
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waves in the lake and imperfections in the glass that covers the 
elevator. 


Figure 3.2 [page 69] Crocuses 
J. Hanan, D. R. Fowler (1990) 


The petals were modeled as Bézier surfaces, with the shapes de- 
termined using Ise. 


Figure 3.4 [page 72] Lily-of-the-valley 
P. Prusinkiewicz, J. Hanan (1987) 


Figure 3.5 [page 74] Development of Capsella bursa-pastoris 
P. Prusinkiewicz, A. Lindenmayer (1987) 


Figure 3.6 [page 75] Apple twig 
P. Prusinkiewicz, D. R. Fowler (1990) 


This twig model was developed in one spring day, looking at a 
real twig nearby. This time is indicative for most inflorescence 
models shown. 


Figure 3.11 [page 81] A mint 
P. Prusinkiewicz (1988) 


Figure 3.14 [page 84] Development of Lychnis coronaria 
P. Prusinkiewicz, J. Hanan (1987) 


Figure 3.17 [page 90] Development of Mycelis muralis 
P. Prusinkiewicz, A. Lindenmayer (1987) 


Figure 3.18 [page 91] A three-dimensional rendering of the Mycelis 
models 
P. Prusinkiewicz, J. Hanan (1987) 


All internodes in the model are assumed to have the same length. 
In reality, the internodes have different lengths, and the structure 
is less crowded. 


Figure 3.19 [page 92] Lilac inflorescences 
P. Prusinkiewicz, J. Hanan, D. R. Fowler 
(1990) 


204 Appendix B. About the figures 


Figure 3.21 [page 94] The Garden of L 
P. Prusinkiewicz, F. D. Fracchia, J. Hanan, 
D. R. Fowler (1988) 


All plants were modeled with L-systems and rendered using the 
IRIS firmware. Images corresponding to different viewing planes 
(the background lilac twigs, the apple twig and the daisies) were 
defocused separately using low-pass filters to simulate the depth 
of field, then composited with a focused image of lilac inflores- 
cences. The sky was generated using a fractal algorithm. 


Figure 3.23 [page 96] Wild carrot 
P. Prusinkiewicz (1988) 


Figure 4.3 [page 102] Close-up of a daisy capitulum 
D. R. Fowler (1988) 


The petals and florets were modeled as Bézier surfaces. 


Figure 4.4 [page 102] Domestic sunflower head 
D. R. Fowler, P. Prusinkiewicz (1989) 


Figure 4.5 [page 105] Sunflower field 
D. R. Fowler, N. Fuller, J. Hanan, A. Snider 
(1990) 


This image contains approximately 400 plants, each with 15 flow- 
ers. A flower has 21 petals and 300 seeds, modeled using 600 tri- 
angles and 400 triangles respectively. Counting leaves and buds, 
the entire scene contains about 800,000,000 triangles. The image 
was ray-traced with adaptive supersampling on a grid of 1024 x 
768 pixels using 45 hours of CPU time on a MIPS M-120 com- 
puter. 


Figure 4.6 [page 106] Zinnias 
D. R. Fowler, P. Prusinkiewicz, J. Hanan, 
N. Fuller (1990) 


The vase was modeled using L.E.G.O. and rendered with a pro- 
cedural texture. The scene was illuminated by one extended light 
source. 
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Figure 4.7 [page 106] Close-up of zinnias 
D. R. Fowler, P. Prusinkiewicz, A. Snider 
(1990) 


This scene was rendered using distributed ray-tracing to simu- 
late the depth field. 


Figure 4.8 [page 108] Water-lily 
D. R. Fowler, J. Hanan (1990) 


Figure 4.9 [page 108] Lily pond 
D. R. Fowler, J. Hanan, P. Prusinkiewicz, 
N. Fuller (1990) 


The wavelets on the water surface were obtained using bump- 
mapping with a procedurally defined texture. 


Figure 4.10 [page 109] Roses 
D. R. Fowler, J. Hanan, P. Prusinkiewicz 
(1990) 


Distributed ray-tracing with one extended light source was used 
to simulate depth of field and create fuzzy shadows. 


Figure 4.11 [page 111] Parastichies on a cylinder 
D. R. Fowler (1990) 


Figure 4.15 [page 116] Pineapples 
D. R. Fowler, A. Snider (1990) 


The image incorporates a physically-based model of a tablecloth 
approximated as an array of masses connected by springs and 
placed in a gravitational field. The scene is illuminated by three 
extended light sources. 


Figure 4.16 [page 117] Spruce cones 
D. R. Fowler, J. Hanan (1990) 


Figure 4.17 [page 117] Carex laevigata 
J. Hanan, P. Prusinkiewicz (1989) 


The entire plant, including the leaves, was modeled using para- 
metric L-systems. 
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Figure 5.2 [page 121] Maraldi figure 
Ned Greene, NYIT (1984) 


The shapes of leaves, calyxes and petals were defined using a paint 
program, by interpreting gray levels as height. Painted textures 
were mapped onto the surfaces of leaves and calyxes. Smooth 
gradation of color across the petals was obtained by assigning 
colors to the vertices of the polygon meshes representing flowers, 
then interpolating colors across polygons using Gouraud shading. 
The vines were rendered with bump-mapping, using a digitized 
image of real bark. 


Figure 5.3 [page 121] The fern 
P. Prusinkiewicz (1986) 


Figure 5.7 [page 125] A rose in a vase 
D. R. Fowler, J. Hanan, P. Prusinkiewicz 
(1990) 


Petals and thorns are Bézier surfaces incorporated into a rose 
model expressed using L-systems. The vase was modeled as a 
surface of revolution. 


Figure 6.3 [page 141] Development of Anabaena catenula 
P. Prusinkiewicz, F. D. Fracchia (1989) 


Each developmental stage is plotted in one scan line. 


Figure 7.13 [page 161] Simulated development of Microsorium 
linguaeforme 
F. D. Fracchia, P. Prusinkiewicz, 
M. J. M. de Boer (1989) 


Cells are represented as polygons, rendered using the IRIS firmware. 
The development can be visualized directly on the screen of an 
IRIS 4D/20 workstation without resorting to single-frame anima- 
tion techniques. 


Figure 7.14 [page 161] Microphotograph of Microsorium linguaeforme 
M. J. M. de Boer, University of Utrecht 


Figure 7.16 [page 163] Simulated development of Dryopteris 
thelypteris 
F. D. Fracchia, P. Prusinkiewicz, 
M. J. M. de Boer (1989) 
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Figure 7.19 [page 169] Developmental sequence of Patella vulgata 
F. D. Fracchia, A. Lindenmayer, 
M. J. M. de Boer (1989) 


Cells are represented as spheres. Intersections of spheres inside 
the modeled embryo are ignored, since they do not affect the 
ray-traced images. 


Figure 7.20 [page 169] An electron microscope image of Patella 
vulgata 
W. J. Dictus, University of Utrecht 


Figure 8.4 [page 180] Fern dune 
P. Prusinkiewicz, D. Hepting (1989) 


The shape of the leaf has been captured using a controlled iter- 
ated function system. A continuous escape-time function defines 
point altitudes, resulting in a surrealistic incorporation of a leaf 
into the landscape. 


Figure 8.9 [page 186] Carrot leaf 
D. Hepting, P. Prusinkiewicz (1989) 


The leaf shape has been modeled using a controlled iterated func- 
tion system. The scene consists of a set of spheres, with the radius 
equal to the distance to the leaf. The image was rendered using 
ray-tracing. 


Figure E.1 [page 191] Water-lilies 
D. R. Fowler, J. Hanan, P. Prusinkiewicz, 
N. Fuller (1990) 


A scene inspired by Water-lilies pool - Harmony in green by 
Claude Monet (1899). All trees and water-lilies were modeled us- 
ing L-systems. The willow twigs bend downwards due to a strong 
tropism effect, simulating gravity. The bridge was modeled using 
L.E.G.O. The sky is a sphere with a procedural texture. The 
entire scene was ray-traced, then the resulting image was repre- 
sented as a set of small circles, with the colors close but not equal 
to the average of pixel colors underneath. This last operation was 
aimed at creating the appearance of an impressionistic painting. 


Figure A.4 [page 200] Virtual lab 
L. Mercer, D. R. Fowler (1990) 


Turtle interpretation of 


symbols 

Symbol Interpretation Page 
F Move forward and draw a line. 7, 46 
i Move forward without drawing a line. 7, 46 
+ Turn left. 7, 19, 46 
— Turn right. Calg 
A Pitch up. 19, 46 
& Pitch down. 19, 46 
\ Roll left. 19, 46 
/ Roll right. 19, 46 
| Turn around. 19, 46 
$ Rotate the turtle to vertical. 57 
[ Start a branch. 24 
| Complete a branch. 24 
{ Start a polygon. 120,127 
G Move forward and draw a line. Do not record a vertex. 122 
Record a vertex in the current polygon. 122. 12% 

} Complete a polygon. 120, 127 
~ Incorporate a predefined surface. 119 
! Decrement the diameter of segments. 26, 57 
: Increment the current color index. 26 


% Cut off the remainder of the branch. 73 
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